Alexis Paljic
email: alexis.paljic@inria.fr

Sabine Coquillart
email: sabine.coquillart@inria.fr

A passive stringed haptic system for immersive environments

We propose a passive haptic feedback system that can provide the user with grounded forces in a 3D manipulation space. A visual/haptic interaction is used to control forces returned by the system, which is composed of cables and a braking system. An experiment, consisting of two sessions of stickiness comparison tasks, is conducted to evaluate the device. The first session is done with the proposed haptic feedback, the second is purely visual, without haptic cues. We propose indicators for task performance and we measure the forces fed back by the device during the experiment. Results show that the device provides the users with actual forces that are consistent with simulated physical properties, and that it enhances user performance compared to the purely visual situation. Finally, the system is installed on a Responsive Workbench on which informal tests show that it provides the users with haptic sensations.

Introduction

A definition of passive haptic systems is systems that are not computer controlled. "Props" [START_REF] Hinckley | Passive real-world interface props for neurosurgical visualization[END_REF] are an example of passive haptic systems; they are rigid or deformable objects which provide feedback simply by their shape or texture, unlike active haptic devices that use computer controlled actuators. Another possible definition was proposed by Swanson [START_REF] Swanson | Implementation of Arbitrary Path Constraints using Dissipative Passive Haptic Displays[END_REF]: "passive haptic displays are systems that are energetically passive, in the sense that they don't add kinetic energy to the system" -unlike active systems-"but are only capable of removing, storing, or redirecting kinetic energy". It is this definition that we will use in this work. Props seem well suited to simulate reaction forces from solid objects such as a virtual cockpit, command buttons or walls, and might also be suitable for compliant objects (sponge, spring). However, it is difficult to change the simulated resistance or compliance of virtual objects since their haptic representation -the physical object-would have to be changed. Some passive devices can adapt to changes in the virtual environment. The PTER (Passive Trajectory Enhancing Robot) was developed by Book et al. [START_REF] Book | The concept and implementation of a passive trajectory enhancing robot[END_REF] and makes use of electromagnetic brakes to constrain user movements. Such systems, that feed back forces to the user by constraining his movements (dissipating his kinetic energy), with computer controlled brakes or clutches are called dissipative passive devices by Swanson [START_REF] Swanson | Implementation of Arbitrary Path Constraints using Dissipative Passive Haptic Displays[END_REF]. Cobots [START_REF] Colgate | Cobots: Robots for collaboration with human operators[END_REF] are steerable systems that can constrain a user's movements by reducing his degrees of freedom. These systems use computer driven actuators, and can be costly and complex to maintain. Pseudo-haptic feedback was proposed by Lecuyer et al. [START_REF] Lecuyer | Pseudo-haptic feedback: Can isometric input devices simulate force feedback?[END_REF], who defined pseudohaptic systems as "systems providing haptic information generated, augmented or modified, by the influence of another sensory modality". The experiment conducted by the authors showed the possibility of providing the user with haptic information by means of a passive input device combined with visual feedback. This passive technique makes no use of computer driven actuators and can be used to simulate different values of mechanical properties such as stiffness. Haptic systems can provide grounded or non-grounded forces. Grounded systems can constrain user movements relative to the ground, such as the Cyber-Force [6], or the Phantom TM [START_REF] Massie | The phantom haptic interface: A device for probing virtual objects[END_REF] and return grounded forces. They are opposed to nongrounded forces that are provided by portable systems, and can constrain the user's movements relative to his own body. The objective of this work is to propose a passive feedback technique that cumulates the advantages of totally passive devices (low cost, no computer driven actuators) and the advantage of computer driven passive devices (adaptation to changes in the virtual environment), that can simulate grounded forces in a 3D workspace. The main question is how to constrain user movements in 3D space, depending on variable forces in the virtual environment, without an active, computer controlled device? The proposal of this work is to use a visual/haptic interaction, inspired by pseudo-haptic feedback. The visual/haptic interaction is meant to induce actions from the user to trigger reaction forces from a passive haptic device based on grounded cables and a braking system. In Section 2, we describe the proposed passive haptic system and evaluate it in Section 3 in a virtual stickiness discrimination task. In Section 4 we discuss the results of the experiment on the basis of user performance and force measurements on the device.

Proposed passive haptic device

Pseudo-haptic feedback, proposed by Lecuyer et al. [START_REF] Lecuyer | Pseudo-haptic feedback: Can isometric input devices simulate force feedback?[END_REF] makes use of an interaction between two sensory modalities, visual and haptic, without using a computer controlled actuator. The authors conducted an experiment using an isometric device (a Spaceball) allowing the user to move a virtual cube on a plane where a "sticky" zone is represented. When the cube is in the sticky zone, the displacement gain (object movement/user movement ratio) of the object is lowered. The reaction of the users was to press the device harder, causing higher reaction force from the device. The authors also conducted an experiment showing that pseudo-haptic feedback can simulate stiffness of springs. For both cases, the device was put on a table and was not meant to move.

The system that we propose is based on a similar interaction between visual and haptic information. Our purpose is to provide the user with grounded forces that can occur in 3D manipulation space in an immersive environment. Thus, the device has to be mobile.

In the experiments on pseudo-haptic feedback, the displacement of the virtual object is rate-controlled, and the reaction of the device to user actions is inherent to its internal stiffness. In the case studied here, the device has to be mobile, and the user control on virtual objects is a position-control. For the rate-control case, control actions cause reaction forces from the isometric device. For the position-control case, there are no such reaction forces. An additional action to the position-control movements of the user is required in order to yield reaction forces from the system. We propose that this additional action of the user be a pressure on the device.

The principle of the system is (see also Figure 3): The user can move a virtual object in 3D space by moving the device. The device contains a tracker for movement acquisition, a force sensor1 and a braking device (see Figure 1). Strings go through the braking element, and when the user presses it, friction forces constrain his movements. When there are no virtual constraints on the object, the user can freely move it by moving his hand (see Figure 3). When the object is under a constraint, the displacement gain (object movement/user movement ratio) is lowered and the object visually moves slowly. In this situation, the user has to press on the device to raise the displacement ratio (according to his pressure measured by the sensor). Consequently, the pressure that the user exerts on the system constrains his translational movements on wires which are attached to the ground. In that situation, he is likely to feel actual resistance forces while he moves in the constrained area, see Figure 3.

In the system that we propose, visual feedback induces the user to press the device depending on visual feedback, and the system returns forces in reaction. We want the user to press on a brake that constraints his own movements. A real life situation where we create ourselves the constraints on our own movements is when we are riding a bicycle: when one presses the brake handles, he actually feels deceleration forces, and is able to modulate the pressure he exerts to brake properly. A one dimensional version of the braking system is shown on Figure 1) and extended to 2D and 3D (Figure 2). This system is based on a visual/haptic interaction: when the object slows, the subject presses the device. Our concern is to see if this passive system actually feeds back forces. In the next Section we describe an experiment meant to investigate this question. When the virtual object (small cube) is under no constraint the user doesn't have to press the braking device to move it, and his movements are free. Right: Constrained displacement. In order to move the virtual object in the sticky area, the user has to press the device. Consequently, his movements along the strings are constrained.

System Evaluation

Method

To evaluate the system, we had a group of subjects perform a stickiness discrimination task. The users had to compare the virtual stickiness of two areas placed on a virtual plane. For this experiment, the simulation was displayed on a 25 inch CRT monitor, and runs on a Sgi Onyx2 computer. Users controlled a virtual cube which could move on the plane and inside the sticky areas. The movements of the cube were visually constrained inside the sticky areas, according to each area's simulated stickiness. The experiment consisted of 2 sessions, one session using the passive haptic feedback (HAPTIC session), and a session with haptic feedback disabled (VISUAL session). The stickiness of the reference area was constant. For the comparison area it varied from -30%, -20%, -10%, 10%, 20% and 30% of the reference stickiness, yielding 6 comparison pairs. Each comparison pair was presented 3 times in a random order, yielding 18 comparison tests for each session. For each test the subject was asked to choose which of the two areas seemed the stickier. 18 subjects participated in the experiment, 16 men and 2 women, aged from 20 to 29. 2 were left handed, 16 right handed.

Gathered data and results

Indicators used to describe user performance were: Score: The overall score is the number of correct answers given by all subjects, considering that the reference for each answer is the relation between the reference and comparison stickiness. Overall score is 69.7% for visual session and 77.4% for haptic session. Just Noticeable Difference : We followed the usual method described in [START_REF] Gescheider | Psychophysics Method, Theory, and Application[END_REF] to compute the Just Noticeable Difference (JND), which is the smallest increment in stimulus (here, comparison stickiness) that is perceived by the user. The PSE is the value of the comparison stimulus that is perceived to be equal to the reference one. JND found is 25.8% for VISUAL session and 22.8% for HAPTIC session. Respective PSE values are 1.13 and 1.23 (virtual stickiness values), corresponding to overestimations of 0.4% and 9% of the reference value (1.125). Applied force : The forces applied on the system were gathered to investigate if actual forces felt by the user are consistent with the simulated stickiness.

Discussion

The main issue in this work was to investigate if the proposed system could return forces that are consistent with the virtual forces we intend to simulate. The data for average applied forces on the comparison area show that the stickier the comparison area, the harder users press the device. Furthermore, this relation between the force we intend to simulate and the average force returned by the system is close to linear (R 2 = 0.9624). Secondly, the system increases performance: this result tends to show that the system brings additional information compared to the visual session. Results for score and applied forces strengthen the idea that this additional information is not an illusion, but indeed actual forces, consistent with the simulated property. We conclude that the system can be used to feed back haptic information: for example it can be suitable to simulate variable electromagnetic force fields in a scientific visualisation application.

Conclusion

We have proposed a passive haptic display that doesn't make use of computer driven actuators and provides the user with grounded forces in 3D workspace. Results show better resolution (JND) in perception of simulated forces when using the system, compared to a situation where no haptic information is available. Secondly, the forces applied by the users to the device are consistent with the simulated property: the stickier a presented virtual area is, the harder the user presses the device and constrains his movements in that area. The system was installed on a Responsive Workbench, showing its adaptability. On the basis of informal tests on the workbench, users reported having haptic sensations such as stickiness or viscosity. This experiment involved comparison between two virtual properties and further work will evaluate the ability of the system to simulate properties that could be compared to real physical properties, for example, present the user with a real viscous object and compare it with a virtual one, simulated with this passive feedback.

Fig. 1 .

 1 Fig. 1. 1D braking solution. When the user applies a pressure p on the braking system, a friction force F occurs during his displacement d. Since the string is fixed to the ground, the friction forces are grounded.

Fig. 2 .

 2 Fig. 2. Left: 2D solution, one end of each string is fixed to the ground, the other is fixed to a winding so the length of the string can vary. Right: Schematic view of the passive haptic system (3D version). 4 strings are attached to 4 corners of a cubic structure, each string can wind and unwind at the opposite corner

Fig. 3 .

 3 Fig.3. Left: Unconstrained displacement. When the virtual object (small cube) is under no constraint the user doesn't have to press the braking device to move it, and his movements are free. Right: Constrained displacement. In order to move the virtual object in the sticky area, the user has to press the device. Consequently, his movements along the strings are constrained.

Proceedings of EuroHaptics 2004, Munich Germany, June5-7, 2004.

(Flexiforce TM ELF TM system, Sensor B101-L [8])


Victoire Marilys

Razakamanana

Martine Audibert
email: martine.audibert@uca.fr

Aina Harimanana

Catherine Araujo Bonjean

Keywords: Madagascar, Malaria, Pneumonia, Community health workers, Difference in difference methodology, Cost-effectiveness I15, I12, I18

In Madagascar, in February 2014, the Ministry of Health and UNICEF implemented a program integrating the diagnosis and treatment of pneumonia into malaria community case management. The objectives of this program were to improve the management of cases of malaria and pneumonia by community health workers to alleviate the problem of accessibility to care and to reduce the number of severe cases treated at health facilities. This paper aims to assess the effectiveness and the efficiency of this. Two districts were taken into account: Andapa received only basic activities ensuring the functionality of the community sites (control district) and Antalaha where all activities related to the program were implemented (treated district). To assess the impact of the program, we use the difference in difference methodology and we compare the period before the implementation of the program in January 2014 and the period when the program is implemented in February 2016. Then cost-effectiveness analysis was made. In Antalaha, although the program has no significant impact on pneumonia, the situation is better than that of Andapa, as in the case of malaria management, the difference is significant between the two districts. The cost-effectiveness analysis also demonstrated that the cost per case of additional malaria and pneumonia treated by CHWs is 2.52 USD (2.44-3.50). However, skills of CHWs should be strengthened especially concerning pneumonia cases management.

Introduction

Since 2000, in Madagascar, significant progress in the reduction of infant mortality was observed. Indeed, between 2000 and 2013, the death rate of children under 5 decreased from 109.2 to 53.4 per 1,000 live births (WDI, 2015). However, this rate remains high compared to the world average of 45.6 in 2013. This rate is 14.3 in Mauritius, a neighboring island; it is 4.4 in France and the lowest infant mortality rate is registered in Norway with 2.8 per 1,000 live births (WDI 2015).

In Madagascar, malaria, pneumonia and diarrhea are the leading causes of death of children. 34% of deaths among children aged between 1 to 59 months are caused by these three diseases (WHO, 2013). However, many deaths could have been avoided through preventive actions, curative treatment and appropriate treatment. The problem is that many households have difficulty access to care because of the geographical distance between their houses and the health center, the shortage of medical staff and the lack of money. Therefore, the: Integrated Community Case Management (ICCM) approach was implemented to solve this problem. Community health workers (CHWs) diagnose and treat simple childhood illnesses. To this end, they receive basic training followed by supervision, tools and inputs for their activities. Thus, at the community level, simple malaria cases are diagnosed by CHWs using Rapid Diagnostic Test (RDT). Positive cases are treated with "Artemisinin-based Combination Therapy" (ACT), simple pneumonia cases are treated with Cotrimoxazole and simple diarrhea cases, with a combination of Zinc-ORS (Oral Rehydration Solution). Severe cases are referred to health centers. According to the World Health Organization (WHO, 2012), this strategy reduces the mortality rate of children under-five and hence induces WHO and UNICEF to promote it [START_REF] Lawn | Donnons sa chance à chaque nouveau-né de l'Afrique : Données pratiques, soutien sur le plan politiques et programmes pour les soins du nouveauné en Afrique, Cape Town : Partnership for Maternal, Newborn and Child Health[END_REF].

In Madagascar, ICCM has been in place in 2011 in all the 22 regions and the CHWs were trained in the management of malaria, pneumonia and diarrhea cases. However, in 2012, inputs stock shortages for the treatment of diarrhea and pneumonia were observed at the community level (PNLP Annual Report, 2012). As acute respiratory infections (ARI), including pneumonia, are the leading causes of death in the Island (MOH, Health Statistics, 2012) a strengthening of the ICCM have been done.

The Malagasy Ministry of Health (MOH) and UNICEF in February 2014 decided to provide to each CHW an ARI-timer ("Acute respiratory infections timer") to detect pneumonia cases and Amoxicillin-DT (dispersible) tablets for the treatment of this disease. The program targets children under-five who are vulnerable to these diseases because of their low immunity (UNICEF, 2014). CHWs' skills were strengthened through training, formative follow-up and retraining. These activities were carried out in Antalaha, a district in the SAVA region. In Andapa, CHWs received initial training and inputs, activities that ensure the survival of sites. Indeed, in areas that have not supported by technical and financial partners (TFP), most community sites are no longer functional. This is the case of the district of Sambava (District Health Service of Sambava, 2015). After the initial training on the implementation of the ICCM in 2011, CHWs in this district have received neither training nor inputs. As a result, they no longer managed cases. Thus, in this district, CHWs do not report their activities to the heads of the health facilities. This situation explains the lack of information and data on the number of cases handled by the CHWs in this district. This region was chosen by the MOH and UNICEF because in 2012, severe malaria is the second leading cause of hospital morbidity in children under-five and pneumonia the fourth (MOH, Health Statistics, 2012). Both diseases are the leading causes of death in this age group (Regional Health Office of SAVA, 2015).

The integration at the community level of pneumonia and malaria management is a new and debate strategy. It implies a reorganization of children healthcare even it is considered as an innovative strategy. In Madagascar, this strategy needs to be boosted as it seems to be low used since its implementation in 2011. UNICEF decided to strengthen this strategy by supporting CHWs training and introducing ARI-timers. But, before embarking on the scaling, evaluate the impact of additional activities, both from the point of view of its effectiveness and its efficiency is needed. It is the objective of our paper.

Literature review

To solve the problem of accessibility to healthcare, CHWs are trained on community case management. According to WHO (2012), this system can reduce costs and decrease the mortality and morbidity rates, especially among children under 5 years old. The study of Theodoratou et al.

(2010) about the case of low income countries shows that malaria community case management can decrease from 60% mortality related to this disease and all-cause mortality in children less than 5 years to 40%. Pneumonia community case management in turn can reduce the mortality related to it by 70%.

In case of fever, CHWs use RDT. A positive test indicates a case of malaria. Simple malaria cases are treated immediately with ACT and severe malaria cases are referred to health facilities. Then, the CHWs count the child's breathing rate with ARI-timer. If the frequency is higher than normal, it is pneumonia. Counting the breathing rate is a problem for some CHWs who have difficulty to assess it. However, drug overdose can cause antibiotic resistance and increase case fatality rate [START_REF] Druetz | The community case management of pneumonia in Africa: a review of the evidence[END_REF]. Some CHWs, specifically the unpaid CHWs, busy with their own daily activities, spend little time on community case management. In Rwanda, for example, they spend 0.2 hour per week (1% of their time) to community case management; the rest is devoted to economic and domestic activities (Jarrah et al., 2013a). In some countries, such as Malawi, CHWs receive wages, so this rate may be as much as 40% of their time [START_REF] Jarrah | Costing of Integrated Community Case Management in Malawi, Analysis report USAID[END_REF].

Like other program, the ICCM can entail direct and indirect costs. These costs vary considerably from country to country. In most cases, the most important expenditures are related to drug supply, training and CHWs equipment (Jarrah et al., 2013a). Table 1 shows the costs of that program in several African countries. They considered two control and two intervention areas. Authors found that the chance not to have a fever in the intervention area is 23.5%. However, this study focused especially on the impact of the use of nets while the program is also related to the training of CHWs and to the improvement of accessibility to ACT. [START_REF] Mubiru | Evaluation of integrated community case management in eight districts of central Uganda[END_REF] assessed the impact of ICCM program in Uganda. From 2010 to 2012 a pilot project has been implemented in eight districts while three other districts were considered as control districts. Mortality rates in the two groups (control and treated) were compared. The method of DID has been made. The impact on mortality has been identified from the method Lives Saved Tool (LIST) 1 . They found that through access to ACT and Amoxicillin, the under-five mortality rate in the intervention area decreased from 50 to 49 deaths per 1,000 live births. However, in the control zone, this rate increased from 63 to 69 deaths per 1,000 live births.

Concerning the assessment of the impact of the ICCM program before 2014, Amouzou et al. ( 2014) identified 11 studies. But they considered only 8 because of the poor quality of the 3 others. All these studies used the DID methodology. Among these studies, only one study has shown that ICCM can reduce infant mortality.

Methodology

Our study was conducted in two (Andapa and Antalaha) of the four districts of the SAVA region, situated in the Northeast of Madagascar.

ICCM Program

The program has been set up in February 2014. All program activities were implemented in Antalaha, chosen by MOH and UNICEF as pilot district. In Andapa (control district), only the basic activities, ensuring the survival of the community sites, were carried out, whose oversight of the drug supply chain. .

The target population is children under-five, residing more than 5 km from the nearest health facility.

Basic activities

In February 2014, cascade training was organized. The officials of the Regional Health Office (RHO) and of the District Health Office (DHS) who had been trained by the MOH officials trained the health workers of health facilities from the districts of Antalaha and Andapa, who, in return, 1 It is a software used to estimate the impacts of maternal and child programs. The user can insert into the tool different information such as demographic details about the country, causes of death by age, the rate of growth, prevalence of Plasmodium falciparum, etc. Then, it builds scenarios. For each scenario, the model produces estimates of deaths avoided and cost estimates.

formed the CHWs on managing malaria and pneumonia cases at community level. Therefore,

CHWs training has been preceded by two orientation meetings and team building. The total number of CHWs trained was: 306 in Andapa and 396 in Antalaha.

Then supervisions of the community pharmacies were carried out regularly. The aim is to implement a system to monitor and control the management of drugs in the region to avoid drug stock-outs in the two districts. These supervisions are also used to check whether the revolving mechanism of drugs is respected and to verify the use of funds generated from the sale.

Program activities in Antalaha

Program activities were to provide formative follow-ups, retrain the CHWs; provide tools such as sheets to note all the cases, sheets for the report as well as start-up kits with ARI-timer and two boxes of Amoxicillin-DT 250 mg; strengthen the community mobilization and home visits. In Andapa (control district), only the basic activities, ensuring the survival of the community sites, were carried out, whose oversight of the drug supply chain.

Refresher training is performed after an evaluation of the level of the CHWs' skills compared to the initial training given to them. Dialogue sessions through home visits where CHWs persuade parents about the need to have their children consulted in the event of illness. CHWs are equipped with consulting cards composed of images on the front and messages on the back. Each CHW has ten days to do the home visits.

Impact analysis: Difference in difference methodology (DID)

In this study, we used DID methodology. The number of cases treated is compared in the following two situations: in which the CHWs benefited from all of the program activities, case of Antalaha (treated group) and in which they received only the core activities ensuring the functionality of their sites, Andapa's case (control group). In order to avoid confusion biases, the introduction of covariates is necessary [START_REF] Gertler | L'évaluation d'impact en pratique[END_REF].

The two strategies we compare are therefore the case of Andapa and that of Antalaha. The aim is to determine whether the additional activities have had an impact on the number of cases of malaria and pneumonia treated by CHWs.

Number of observations

Normally, the number of observations and the number of health facilities are the same since CHWs reports are compiled at the health facilities. However, only 15 out of the 28 health facilities of Andapa and 29 out of the 32 Antalaha health facilities provided reports. The control group and the treated group are therefore composed respectively of 15 and 29 health facilities. The others are mainly those suffering from problems of geographical accessibility. For example two health facilities of Antalaha are at four days walk from the last place accessible by motorcycle. The heads of these health facilities rarely come to the DHS.

Model and variables description

We use the method of [START_REF] Card | Minimum Wages and Employment: A Case Study of the Fast-Food Industry in New Jersey and Pennsylvania[END_REF]. Their study is about the impact of rising wages on the level of employment in fast food in New Jersey. To do this, they compare the situation in New Jersey (treated group) with that in Pennsylvania where there was no wage increase (control group).

Here, we compare the number of cases of malaria and pneumonia treated by the CHW in Andapa district with those in Antalaha district.

Let:

[X] is the access to the program. The value of this random variable is 1 for Antalaha and 0 for Andapa.

[Y] represents the effectiveness of the program. It expresses the number of cases of malaria and pneumonia treated by CHWs divided by the number of the target population. This indicator is the equivalent of the consultation rate in the case where the caregiver is a CHW. It should be higher in the case of the treated district than in the case of the control district.

[Y 0 ] is the number of cases handled by CHWs before the program implementation in t 0 . [Y 1 ] is the number of cases treated by CHWs in t 1 , period when the population appropriates the program. ∆Y represents the change in the number of cases that CHWs have managed between t 1 and t 0 .

Tr is the treated group and C the control group. DID estimation is:

) ( ) ( C it Tr it Y E Y E     
, with E the effectiveness.

The average effect of the program in Antalaha is:

) 1 ( 0 1     X Y Y E Y TT
To avoid the risk of selection bias, the situation before the program implementation must be the same for the two districts.

Then, the link between results and participation in the program could be explained by other variables. The introduction of covariates can also influence outcomes [START_REF] Khandker | The International Bank for Reconstruction and Development[END_REF]. This approach can avoid confusion bias. Therefore, the number of CHWs attached to each health facility [CHW], the distance between each community site and the district to which it is attached [KM], as well as the prevalence of malaria and pneumonia [Prevalence] at the level of each municipality of the health facilities are considered.

Source of data about effectiveness

Health data came from the District Health Service (DHS) where the monthly activity reports of the health facilities and those of CHWs are collected. These reports contain information on the number of children in care per month by each health facility and by CHWs. They include, among others the number of fever cases, cases tested by RDT, malaria and pneumonia cases and the number of children referred to health facilities. Monthly meetings are held at health facilities during which CHWs must submit activity reports concerning the health situation at their respective sites.

Health facilities should have all the data concerning CHWs activities. However, some CHWs are not attending the meetings mainly because of the problems of accessibility. Indeed, transportation is not covered by the program, whereas in order to attend meetings, some CHWs have to pay travel fare which is equivalent of 3.5 USD (survey, 2016). Note that for those who practice in very remote sites, carts are the only available means of transportation, so, data are incomplete.

Cost estimation

As ICCM induces a reorganization of child health care, and may be some constraint for health staff, we decided to include both direct (financial) and indirect (economic) costs. The total cost is the total expenditure in the districts of Antalaha and Andapa for the period from February 2014, the date of the program implementation, to February 2016, period after the implementation.

Efficiency: Cost-effectiveness analysis

Program efficiency is assessed using cost-effectiveness analysis. This is done by comparing the strategies and evaluating the incremental cost for an additional effectiveness (O'Brien et al., 1994).

In the studied case, the aim is to verify whether carrying out additional activities in addition to the core activities is expensive or not in comparison to the obtained results. To this end, the case of Andapa is compared, where only the basic activities relating to community health have been carried

out with that of Antalaha where all the activities related to the program have been implemented.

According to Drummond, Sculpher and Torrance (2005), the reference is the incremental costeffectiveness ratio (ICER) where the numerator is the difference of the costs between the two strategies, and the denominator the difference between their effectiveness.

Results

Descriptive analysis: Evolution of the number of children managed by

CHWs and treated at health facilities in Antalaha and Andapa

In the Figure 1, the x axis is the period (month) and the y-axis is the number of cases divided by the number of target children. It should be noted that the average number of episodes of malaria and pneumonia is 1 episode per child per year in Antalaha and Andapa (Survey, 2016).

Thanks to immediate care at community level, severe cases received by health facilities should be reduced. Table 2 shows the variation in the number of pneumonia and malaria cases treated at health facilities. In February 2016, the number of malaria cases treated at health facilities in the two districts decreased. This decline is greater for the case of Antalaha than that of Andapa.

Concerning pneumonia, the number of consultations decreased in Andapa but significantly increased in the health facilities of Antalaha.

Difference between Andapa and Antalaha effectiveness

In the following table (Table 3), model (I) is the results of a simple DID. Model (II) is obtained after the introduction of covariates. In model (III), DID is associated with propensity score matching. The propensity score is the probability for an individual i to receive treatment considering his characteristics V i [START_REF] Khandker | The International Bank for Reconstruction and Development[END_REF]. The logit model is as follows:

) 1 ( ) ( i i i i v V Z P x e   
. The aim is to correct the heterogeneity between the treated and the control groups and to avoid an overestimation of the impact of the program. Source: Authors, from DHS data A significant difference was observed between the number of malaria cases managed by CHWs in Andapa and Antalaha between t 1 and t 0 . But the difference is not significant for the case of pneumonia.

Program costs

The total cost of interventions is twice higher in Antalaha than in Andapa (Table 4). Furthermore, the cost of the program per target child for each of the two scenarios (Andapa and Antalaha) belongs to the cost range estimated by [START_REF] Collins | The cost of integrated community case management (iCCM) programs: a multi-country analysis[END_REF] studying the cost of the same program in seven sub-Saharan Africa countries (Table 1). For the estimation of efficiency previously, we use the results of DID model (III), when we use those of model (II), the ICER is 3.50 USD and 2.44 USD with the model (I). Therefore, the cost range belongs to (2.44 USD -3.50 USD).

Discussion

The results of the DID estimation show that the difference of the number of cases treated by CHWs in Andapa and Antalaha during t 1 and t 0 is significant for malaria cases and not for pneumonia cases. The increase of the number of pneumonia consultations in Antalaha health facilities can be explained by the fact that CHWs are not adequately trained on the prevention and control of this disease and they do not have the required skills to do this. Indeed, according to the report on formative follow-up of September 2015, 23% of CHWs faced a problem to count respiratory rate and 16% made errors in treatment (DHS Antalaha, 2015).

One of the problems with the program is that some CHWs give Amoxicillin-DT to patient with cough or simple cold (Survey, 2016). Such practices increase resistance to antibiotics and make the use of health facilities inevitable. However, concerning malaria management, the difference between the number of cases treated by CHWs in Andapa and Antalaha is significant. Therefore, Andapa should receive the same activities as Antalaha. But the program must be reviewed. The incremental cost for an additional efficiency is 2.52 USD: a cost which is not very high. So to ensure the quality of the community case management, additional activities can be implemented. For example, we can increase the number of retraining.

Thanks to the program, household health expenditures can decrease. The program implementation is the responsibility of the MOH or the TFP, so, it is necessary to know their willingness to pay.

However, it was difficult to obtain this information. Indeed, the share of the government budget for health is already insufficient (MOH, 2014) (8.06% of the total budget in 2012, 7.10% 1 in 2013 and 9.62% in 2014).

Conclusion

In February 2014, MOH and UNICEF implemented a program integrating the diagnosis and treatment of pneumonia into malaria community case management. The objectives of this program were to improve the management of cases of malaria and pneumonia by CHWs to alleviate the problem of accessibility to care and to reduce the number of severe cases treated at health facilities.

In this case, CHWs should be able to manage adequately simple malaria and pneumonia cases to win local population trust. To this aim, two districts of the SAVA region were considered: Andapa received only basic activities ensuring the functionality of the community sites and Antalaha where all activities related to the program were implemented. For the analysis, Andapa was taken as the control district and Antalaha as the treated district.

The DID estimation shows that the difference of the number of malaria cases treated by CHWs in Andapa and in Antalaha is significant before and after the program. However, concerning pneumonia cases treated by CHWs, this difference is not significant. Then, it was found that in February 2016, the number of malaria cases treated at health facilities decreased and this decrease was more significant in Antalaha than in Andapa. However, concerning pneumonia, the number of cases treated at health facilities decreased in Andapa and increased significantly in Antalaha.

Therefore, skills of CHWs should be strengthened especially concerning pneumonia cases management.

In Antalaha, although the program has no significant impact on pneumonia, the situation is better than that of Andapa, as in the case of malaria management, the difference is significant between the two districts. Our results showed that the costs of this program are within the cost range estimated by [START_REF] Collins | The cost of integrated community case management (iCCM) programs: a multi-country analysis[END_REF] studying the cost of the same program in seven sub-Saharan Africa countries. The cost-effectiveness analysis also demonstrated that the cost per case of additional malaria and pneumonia treated by CHWs is 2.52 USD (2.44-3.50). However, it has been difficult to define a threshold beyond which this intervention will not be cost-effective. But as said by WHO (2003), the priority is to ensure equity and to ensure that everyone has access to care.

Financial costs are those

  related to training, supervision, upgrading, supply and distribution of inputs and outputs in the DHS reports. Also input costs: ARI-timer, medicines and management tools costs are included. Indirect costs are those related to time. Data come from surveys of each individual responsible of the program. Thus, the technical and financial partners, officials in the MOH, at the RHO, at the DHS of Antalaha and Andapa, at the health facilities and CHWs were interviewed to determine the time spent on the implementation of the program. Different questionnaires 2 were used according to the targeted officials. It is obtained by multiplying the time devoted to the program by the hourly wage.

  costs in Andapa C2: program costs in Antalaha E1: intervention effectiveness in Andapa E2: intervention effeciteveness in Antalaha In our case, the effectiveness is the number of additional cases treated by CHWs thank to the new strategy. It is obtained by comparing the number of cases treated by CHWs before the implementation of the program t 0 in January 2014 with that registered in t 1 in February 2016. We did not use a threshold, because of the difficulty to have information about willingness to pay.

Figure 1 :

 1 Figure 1: Evolution of the number of children managed by CHWs in Antalaha and Andapa (divided by the number of target children) between February 2014 and February 2016

Table 1 :

 1 ICCM costs in other countries

	Authors	Countries	Costs	Recommandations
	Jarrah et al. (2013) Rwanda	0.50 USD/total	-Find means to reduce
			population /year	supervision costs
			2.22 USD/target	-Avoid the problems of
			children/ year	registers and drugs stock-
				outs
	Jarrah et al. (2014) Sierra Leone	10.2 USD/ target	-Find means to reduce
			children/year	supervision costs
	Collins et al.	Cameroon, the	2.07 to 10.26 USD/	-Find means to reduce
	(2014)	Democratic	target children/year	supervision costs and to
		Republic of the		ensure the quality of care.
		Congo, Malawi,		-Avoid bottlenecks like
		Senegal, Sierra		drugs stock-outs
		Leone, South		-Involve community
		Sudan and		leaders
		Zambia
	Source: Authors, 2016

Concerning the impacts of the program, Abegunde et al. (2016) have conducted an impact assessment of the ICCM in Nigeria. Difference in difference (DID) methodology and logistic regression were made. The period before and after the program, 2013 and 2015, were compared.

Table 2 :

 2 Evolution of the number of malaria and pneumonia cases treated at health facilities in

	February 2014 (t 0 ) and February 2016 (t 1 )
		Malaria	Pneumonia
		Antalaha Andapa	Antalaha	Andapa
	t0	3.54	1.28	0.13	0.40
	t1	2.23	1.10	0.70	0.21
	Différence t1-t0	-1.31	-0.17	0.57	-0.19
		Source: Authors, from RHO data

Table 3 :

 3 DID estimations

	Outcome		Malaria/target population (I) (II) (III)	Pneumonia/target population (I) (II) (III)
	T0	Control group	0.20	0.10	0.14	0.46	1.23	0.37
		Treated group	0.44	0.25	0.44	0.82	1.82	0.82
		Difference	0.24	0.15	0.30	0.36	0.59	0.45
			(0.26)	(0.46)	(0.11)	(0.39)	(0.15)	(0.22)
	T1	Control group	0.13	-0.02	0.10	0.17	1.08	0.11
		Treated group	0.99	0.63	0.99	1.42	2.23	1.42
		Difference	0.86***	0.65***	0.89***	1.25***	1.15***	1.31***
			(0.00)	(0.00)	(0.00)	(0.00)	(0.00)	(0.00)
	Double difference	0,62**	0.50*	0.59**	0.89	0.56	0.87
			(0,04)	(0.08)	(0.02)	(0.13)	(0.34)	(0.10)
	Observations	88	88	88	88	88	88
		CHW	-	0.00	0.03	-	-0.07***	0.04
	Covariates	KM	-	(0.66) -0.00	(0.38) 0.01*	-	(0.00) -0.00	(0.22) 0.02*
				(0.79)	(0.08)		(0.23)	(0.07)
		Prevalence	-	0.01***	-0.02		0.02**	-0.02
				(0.01)	(0.71)		(0.03)	(0.31)
			*, **, ***: significant respectively at 10%, 5% and 1%

Table 4 :

 4 Total costs and mean costs in USDThe program implemented in Antalaha allows the diagnosis and treatment at Community level 76 830 malaria and pneumonia cases in comparison of the number of cases treated in Andapa for an additional cost of 193,647.35 USD. This is equivalent to 2.52 USD per additional malaria or pneumonia case (Table5).

	Costs		Andapa	Antalaha
	Activities		37,418.00	130,347.97
	Inputs		71,806.84	114,999.52
	Indirect costs	27,447.29	84,971.99
	Total		136,672.13	330,319.48
	Number of trained CHWs	306	396
	Number of target population	39,678	56,458
	Costs/CHW	447.71	1,017.4
	Costs/target population	3.44	5.85
		Source: Authors, 2016
	4.4.	Incremental cost-effectiveness ratio estimation (ICER)
	Table 5: ICER estimation (USD)
			Malaria + pneumonia
		C2 -C1	193,647.35
		E2 -E1	76,830
		ICER	2.52
		Source: Authors, from DHS data

Études et Documents n° 6, CERDI, 2017

Available upon request. Études et Documents n° 6, CERDI, 2017

Acknowledgments

We are grateful to the Regional Health Officials of SAVA and the District Health Service officials of Antalaha and Andapa for the data. We thank Rakotobe V.M.X. for her proofreading. Finally, many thanks to CERDI and Pasteur Institute of Madagascar.

Marilys Victoire Razakamanana

This work was supported by the LABEX IDGM+ (ANR-10-LABX-14-01) within the program "Investissements d'Avenir" operated by the French National Research Agency (ANR). http://www.cerdi.org/


Youcef Mammeri
email: youcef.mammeri@u-picardie.fr

Damien Sellier
email: damien.sellier@scionresearch.com

A surface model of nonlinear, non-steady-state phloem transport

Keywords: Phloem transport, Diffusion-convection equation, Well-posedness, Finite elements. MS Codes. 35K61, 65M60, 92C80

Phloem transport is the process by which carbohydrates produced by photosynthesis in the leaves get distributed in a plant. According to Münch, the osmotically generated hydrostatic phloem pressure is the force driving the long-distance transport of photoassimilates. Following Thompson and Holbrook[1]'s approach, we develop a mathematical model of coupled water-carbohydrate transport. It is first proven that the model presented here is well posed and preserves the positivity. The model is then applied to simulate the flow of phloem sap for an organic tree shape, on a 3D surface and in a channel with orthotropic hydraulic properties. Those features represent an significant advance in modelling the pathway for carbohydrate transport in trees.

Introduction

Expanded knowledge of the carbohydrate pathway in trees is critical in agriculture, forestry and ecology. Improving on our understanding of translocation is particularly important in a climate change context so as to anticipate the effects of future environmental conditions on tree growth and carbon sequestration. Phloem transport is the process by which carbohydrates produced by photosynthesis in the leaves get distributed within a plant. Efficient transport ensures that the carbohydrate requirements of living tissues (respiration, growth) are met throughout the organism. In most trees, the phloem is a tissue layer located under the bark. Phloem is very thin (typically from 0.5 to 5 mm), i.e. at least 2 orders of magnitude less than its other dimensions. Because of that, it can be described as a three-dimensional manifold with a shape closely matching tree's external shape (minus the offset of bark's thickness). In a schematic view, sap flows in the phloem from leaves (sources) to roots (sinks). In reality, sinks are not only located at one end but also distributed all along the pathway. Even leaves can act as carbohydrate sinks when deciduous trees initiate new leaf growth at the beginning of the growing season. Phloem sap mostly consists in water carrying photosynthates by bulk. Inside the phloem, the sap moves through a network of elongated and interconnected cells. Those cells are referred to as sieve tube elements in angiosperms and as sieve cells in conifers. According to Münch [START_REF] Münch | Die Stoffbewegungen in der Pflanze[END_REF], the osmotically generated hydrostatic phloem pressure is the force driving the long-distance transport of photoassimilates. Münch's hypothesis is generally accepted. For over 60 years, it has served as the basis for mathematical models of the phloem transport process [START_REF] Horwitz | Some simplified mathematical treat-ments of translocation in plants[END_REF][START_REF] Eschrich | Solution flow in tubular semipermeable membranes[END_REF][START_REF] Christy | A mathematical treatment of Münch's pressure-flow hypothesis of phloem translocation[END_REF][START_REF] Ferrier | Further theoretical analysis of concentrationpressureflux waves in phloem transport systems[END_REF][START_REF] Lang | A model of mass flow in the phloem[END_REF][START_REF] Goeschl | Physiological implications of the MünchHorwitz theory of phloem transport: effect of loading rates[END_REF][START_REF] Phillips | Asymptotic analysis of flow in sieve tubes with semi-permeable walls[END_REF][START_REF] Thompson | Application of a Single-solute Non-steady-state Phloem Model to the Study of Long-distance Assimilate Transport[END_REF][START_REF] Jensen | Optimality of the Mnch mechanism for translocation of sugars in plants[END_REF][START_REF] Cabrita | Hydrodynamics of steady state phloem transport with radial leakage of solute[END_REF].

We present the mathematical foundations and an implementation for a surface dynamic, anisotropic model of phloem transport. The purpose of this model is not to elaborate on the finer points of the physical process underlying phloem sap flow such as the contribution of diffusion [START_REF] Payvandi | Mathematical modelling of the Phloem: the importance of diffusion on sugar transport at osmotic equilibrium[END_REF], the radial leakage of solutes [START_REF] Cabrita | Hydrodynamics of steady state phloem transport with radial leakage of solute[END_REF] or the presence of a relay mechanism [START_REF] Hölttä | Linking phloem function to structure: analysis with a coupled xylemphloem transport model[END_REF]. Those aspects can be investigated individually. We adopt a standard mathematical model of coupled water-solute transport that is similar to Thomson and Holbrook's [START_REF] Thompson | Application of a Single-solute Non-steady-state Phloem Model to the Study of Long-distance Assimilate Transport[END_REF]. In this study, the model has been developed with large scale domains and arbitrary geometry in mind. In previous modelling attempts, phloem transport is treated as a one-dimensional process with the sap flowing though a file of sieve-tube elements connected end-to-end. Because plant stems are generally elongated with a high aspect ratio, they appear to be tubes or pipes. As a result, the 1D approach is also employed to simulate transport at the scale of the entire organism [START_REF] Hölttä | Modeling xylem and phloem water flows in trees according to cohesion theory and Münch hypothesis[END_REF][START_REF] Lacointe | Modelling phloem and xylem transport within a complex architecture[END_REF]. However, a 1D model implicitly assumes that carbohydrates are in a common pool at any given position along stem's length. In trees, that assumption may not be valid.

Carbohydrates produced by a single branch are translocated along a downward, helicoidal pathway on the stem surface [START_REF] Hansen | The fate and path of assimilation products in the stem of 8-year-old Scots pine (Pinus sylvestris L.) trees[END_REF][START_REF] Kagawa | CO 2 pulse-labelling of photoassimilates reveals carbon allocation within and between tree rings[END_REF]. That singular trajectory highlights two key points: i) carbohydrates predominantly follow the orientation of sieve elements with little lateral dispersion and ii) the direction of translocation does not correspond to the long axis of shoots and roots. In that context, most of the carbohydrates produced by a source are only available to sinks with a direct hydraulic connection to the source. Therefore, the difference in hydraulic properties along and across sieve elements as well as the lateral positioning of sources and sinks are essential to understand phloem transport in trees. Taking those features into account can be achieved by describing transport as a two-dimensional process [START_REF] Sellier | Phloem transport in trees: A generic surface model[END_REF]. The model we present is a surface model in that the flow of water and solute is neglected within the thickness of the phloem (i.e. intra-phloem radial transport). In mathematical terms, it does not present additional difficulty to extend the current model to the full 3D case. However, there are several reasons not to do so. Firstly, the macroscopic approach used here may not be applicable to a tissue less than a dozen cells wide [START_REF] Gričar | Number of cells in xylem, phloem and dormant cambium in silver fir (Abies alba), in trees of different vitality[END_REF]. Secondly, the fact that width is several orders of magnitude less than other conduit dimensions and that the computational domain must be defined for its smaller dimension would make the problem numerically untractable. Thirdly and more importantly, radial transport appears to follows a different cellular pathway, through ray parenchyma [START_REF] Spicer | Symplasmic networks in secondary vascular tissues: parenchyma distribution and activity supporting long-distance transport[END_REF][START_REF] Pfautsch | Phloem as Capacitor: Radial Transfer of Water into Xylem of Tree Stems Occurs via Symplastic Transport in Ray Parenchyma[END_REF], and Münch's hypothesis may not apply. On the other hand, using a surface model is not incompatible with simulating transport for three-dimensional surfaces and describing radial flow to adjacent tissues (as boundary conditions).

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

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

Model description

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

Figure 1: Schematic representation of a tree and the layered organisation of its secondary tissues: phloem, vascular cambium and xylem (bark not shown). Flows of water/solute in a phloem element also shown.

Equations statement

The transport of a volume θ of water in a surface Ω, a bounded domain of R 2 of thickness e, is defined by the mass balance conservation

∂ t θ + ∇ • J θ + H θ = 0,
where J, the water flux, is given by Darcy's law

J θ = - e µ (k∇P ) .
Here P is phloem's pressure in Pa, µ is the sap viscosity in Pa s and depends on the concentration C in solute, and k in m 2 is the orthotropic permeability matrix

k = k x 0 0 k y .
Typically, the C value is 20% (≈ 630 mol m -3 ). At that C value, µ is twice that of pure water.

C values up to 1000 mol m -3 and above are physiologically realistic. The dependency of µ on C also theoretically affects flow efficiency [START_REF] Lang | A model of mass flow in the phloem[END_REF]. It cannot be neglected beforehand.

H θ denotes the radial water flux at the boundary between phloem and xylem. Here, the radial flux is function of the differential of water potential between phloem and xylem [START_REF] Thompson | Application of a Single-solute Non-steady-state Phloem Model to the Study of Long-distance Assimilate Transport[END_REF][START_REF] Hölttä | Modeling xylem and phloem water flows in trees according to cohesion theory and Münch hypothesis[END_REF]:

H θ = -L R (ψ -P + RT C) -V s U,
and L R the radial hydraulic conductivity in m Pa -1 s -1 , ψ the xylem hydrostatic pressure, R the gas constant (in J mol -1 K -1 ), T the temperature in K, V s the partial molal volume of sucrose in m 3 mol -1 and U the sucrose unloading (radial). It translates boundary conditions on entrance and exit as source and sink terms as follows

U =    U in Ω l the loading area (source) -U C * C in Ω u the unloading area (sink) 0 elsewhere.
Here U denotes a constant loading rate (mol m -2 s -1 ) and C * a reference sucrose concentration (mol m -3 ). On the other hand, the variation of the volume θ depends on the pressure via the phloem thickness (e in m) and phloem Young's modulus (E in Pa) as

∂ t θ = e E ∂ t P.
In other terms, the phloem is deformable and thickness depends on P . Concerning the amount of sucrose, the concentration C is governed by

e∂ t C + ∇ • J C + H C = 0,
where H C = U and J C = CJ θ . To sum up, the mass balance conservation is written, for x ∈ Ω and t > 0

e E ∂ t P -∇ • e µ (k∇P ) = L R (ψ -P + RT C) + V s U (1.1)
e∂ t C -∇ • e µ C (k∇P ) = U, (1.2)
with initial data P (0, x) = P 0 (x) and C(0, x) = C 0 (x).

(

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

∂ n P (t, x) = ∂ n C(t, x) = 0.
(1.4)

Theoretical qualitative study

Since the sap viscosity and the phloem thickness can depend on the concentration of sucrose we assume in the following that there exist four positive constants µ 1 , µ 2 , e 1 , e 2 such that 0 < µ 1 < µ(C) < µ 2 and 0 < e 1 < e(P ) < e 2 .

To deal with the classical solution, the main difficulty is the lack of regularizing term in the mixed parabolic-hyperbolic system for the second variable C [START_REF] Kreiss | Initial-Boundary Value Problems and the Navier-Stokes Equations[END_REF]. Nevertheless, the parabolic regularization, for ε > 0, is well posed

∂ t P -κ∆P + ηP = F (∇P, C) (1.5) ∂ t C -ε∆C + γC = G(∇P, C, ∇C), (1.6)
where

F (∇P, C) = ∇ kE µ • ∇P + EL R e (ψ + RT C) + EV s U e G(∇P, C, ∇C) = ∇ • 1 µ C (k∇P ) + U e . Theorem 1.1 Let (P 0 , C 0 ) ∈ W 1,∞ (Ω) × L ∞ (Ω)
be positive initial data. Then there exists T > 0 and a unique weak solution

(P, C) ∈ L ∞ (0, T ; W 1,∞ (Ω)) × L ∞ (0, T ; L ∞ (Ω)) of initial-boundary values problem (1.5)-(1.6).
Proof. From Theorem 6.19 of [START_REF] Ouhabaz | Analysis of Heat Equations on Domains[END_REF], G the kernel of the heat equation with homogeneous Neumann boundary conditions defined by

∂ t G(t, x, y) -κ∆G(t, x, y) + ηG(t, x, y) = 0 ∂ n G(t, x, y) = 0 for x ∈ ∂Ω, and
∀(x, y) ∈ Ω 2 , lim t→0 + G(t,
x, y) = δ y (x) satisfies for all y ∈ Ω, there exists a constant c Ω > 0, depending on Ω, such that for 0 < t < T

||G(t, ., y)|| L 1 (Ω) ≤ c Ω , ||∇ x G(t, ., y)|| L 1 (Ω) ≤ c Ω max(1, t -3/4 ), ||∇ y G(t, ., y)|| L 1 (Ω) ≤ c Ω max(1, t -3/4 ). We deduce for u ∈ L ∞ (0, T ; L ∞ (Ω)) Ω G(t, x, y)u(t, y)dy L ∞ (Ω) ≤ c Ω sup t∈[0,T ] ||u(t, .)|| L ∞ (Ω) t 0 Ω
G(t -s, x, y)∇u(t, y)dy

L ∞ (Ω) = t 0 Ω ∇ y G(t, x, y)u(t, y)dy L ∞ (Ω) ≤ c Ω max(T, T 1/4 ) sup t∈[0,T ] ||u(t, .)|| L ∞ (Ω) .
According to the Duhamel formula, (P, C) is solution of the initial-boundary values problem (1.5)-(1.6) if and only if (P, C) is a fixed point of

Φ(P, C)(t, x) := Ω G(t, x, y)P 0 (y)dy + t 0 Ω G(t -s, x, y)F (s, y)dyds, Ω G ε (t, x, y)C 0 (y)dy + t 0 Ω G ε (t -s, x, y)G(s, y)dyds in the functional space X T = B T,1 × B T,0 where B T,s := u ∈ L ∞ (0, T ; W s,∞ (Ω)); sup t∈[0,T ] ||u(t, .)|| L ∞ (Ω) ≤ M . Let (P, C) ∈ X T . We have, for t ∈ [0, T ] Ω G(t, x, y)P 0 (y)dy + t 0 Ω G(t -s, x, y)F (s, y)dyds L ∞ (Ω) ≤ c Ω ||P 0 || L ∞ (Ω) + c Ω max(T, T 1/4 ) sup t∈[0,T ] ||P (t, .)|| L ∞ (Ω) + T sup t∈[0,T ] ||C(t, .)|| L ∞ (Ω) + T , ∇ x Ω G(t, x, y)P 0 (y)dy + t 0 Ω G(t -s, x, y)F (s, y)dyds L ∞ (Ω) ≤ c Ω ||∇ x P 0 || L ∞ (Ω) + c Ω max(T, T 1/4 ) sup t∈[0,T ] ||∇ x P (t, .)|| L ∞ (Ω) + sup t∈[0,T ] ||C(t, .)|| L ∞ (Ω) + 1
and using Young's inequality

Ω G ε (t, x, y)C 0 (y)dy + t 0 Ω G ε (t -s, x, y)G(s, y)dyds L ∞ (Ω) ≤ c Ω ||C 0 || L ∞ (Ω) + c Ω max(T, T 1/4 ) sup t∈[0,T ] ||C∇P (t, .)|| L ∞ (Ω) ≤ c Ω ||C 0 || L ∞ (Ω) + c Ω 2 max(T, T 1/4 ) sup t∈[0,T ] ||C|| L ∞ (Ω) + ||∇P (t, .)|| L ∞ (Ω) 2 .
In summary,

||φ(P, C)|| X T ≤ c Ω ||P 0 || L ∞ (Ω) + ||∇P 0 || L ∞ (Ω) + ||C 0 || L ∞ (Ω) + c Ω max(T, T 1/4 ) 1 + ||(P, C)|| X T + ||(P, C)|| 2 X T Choosing M = 2c Ω ||P 0 || L ∞ (Ω) + ||∇P 0 || L ∞ (Ω) + ||C 0 || L ∞ (Ω) , it comes ||φ(P, C)|| X T ≤ M 2 + c Ω max(T, T 1/4 ) 1 + M + M 2
and φ maps X T onto itself as soon as max(T,

T 1/4 ) ≤ M 2cΩ(1+M +M 2 )
. Similarly, for (P 1 , C 1 ) and (P 2 , C 2 ) in X T , we obtain

||φ(P 1 , C 1 ) -φ(P 2 , C 2 )|| X T ≤ c Ω max(T, T 1/4 ) (||(P 1 , C 1 ) -(P 2 , C 2 )|| X T + sup t∈[0,T ] ||C 1 ∇P 1 (t, .) -C 2 ∇P 2 (t, .)|| L ∞ (Ω) . Since C 1 ∇P 1 -C 2 ∇P 2 = (C 1 -C 2 )∇P 1 + C 2 ∇(P 1 -P 2 ), it provides ||φ(P 1 , C 1 ) -φ(P 2 , C 2 )|| X T ≤ c Ω M max(T, T 1/4 )||(P 1 , C 1 ) -(P 2 , C 2 )|| X T ,
and φ is a contraction mapping if max(T, T 1/4 ) < 1 cΩM .

The conservation of the positivity is needed.

Proposition 1.2 Let (P 0 , C 0 ) ∈ W 1,∞ (Ω) × L ∞ (Ω) be positive initial data.
Then, for all time t ≥ 0 and a.e. x ∈ Ω, P (t, x) ≥ 0 and C(t, x) ≥ 0.

Proof. We have

∂ t P -∇ • E µ (k∇P ) = f (P, C), with f (P, C) := EL R e (ψ -P + RT C) + EV s U e = (f (P, C) -f (0, C)) + f (0, C) = ∂f ∂P (θP, C) + f (0, C),
where 0 ≤ θ ≤ 1. Consider Q = exp(λt)P , the equation becomes

∂ t Q -∇ • E µ (k∇Q) -λ + ∂f ∂P (θP, C) Q = exp(-λt)f (0, C). (1.7)
If there exists (x * , t * ) such that Q(t * , x * ) = min t,x Q(t, x) < 0, then

∂ t Q(x * , t * ) = 0, ∇Q(x * , t * ) = 0, ∆Q(x * , t * ) ≥ 0,
and λ can be chosen such that the left-hand size of (1.7) and the right-hand size have opposite sign.

We deal similarly for C.

The phloem fills up and empties with respect to in and out flows. If V(Ω l ) = V(Ω u ), then the total amount of sucrose Ω C(t, x)dx is increasing as soon as C < C * in Ω u , whereas it becomes decreasing when C > C * in Ω u . More precisely, we have:

Proposition 1.3 For all time t ≥ 0 d dt Ω C(t, x)dx = U e V(Ω l ) - Ωu C(t, x) C * dx .
Proof. Integrating (1.2) over Ω gives

d dt Ω Cdx = Ω ε∆C + ∇ • 1 µ C (k∇P ) + U e dx = Ω U e dx = U e V(Ω l ) - Ωu C C * dx .
Remark 1.4 Let C * > 0 be a given unloading strength. Then the steady state solution C satisfies

Ωu C(x)dx = C * V(Ω l ).
Remark 1.5 It seems difficult to obtain uniform estimates independent of ε to evaluate the limit as ε → 0. We plan to study this issue in a future work. Nevertheless, the positivity of the pressure P and the mass behavior Ω C remain true when ε = 0.

Algorithm framework

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

Let n ∈ N * . P n , respectively C n , denotes the approximation at time t n = n∆t of the pressure P , respectively the carbohydrate concentration C.

Splitting

To deal with the nonlinearity, the equation (1.2) is successively split [START_REF] Hundsdorfer | Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations[END_REF] with a first order in time as, for x ∈ Ω and t ∈ [t n , t n+1 ]

∂ t C - k µ ∇P • ∇ C - U e = 0 and ∂ t C -C ∇ • k µ ∇P = 0, with C(x, t n ) = C(x, t n+1
). Since C(x, t) > 0 as soon as C 0 (x) > 0, the change of variables C = log(C) is applied to the second equation to obtain

∂ t C -∇ • k µ ∇P = 0.
Finally, the transformation C = exp( C) enforces the positivity.

Space and time discretization

Let ϕ 1 , ϕ 2 , ϕ 3 be test functions, with ϕ i , ∇ϕ i in L 2 , for 1 ≤ i ≤ 3. Finite elements method in space is used, while the time discretization is obtained with a semi-implicit scheme.

Algorithm 1 Semi-implicit scheme Given P 0 , C 0 and N ∈ N * .

for n = 1 to N do

C n+1/2 -C n ∆t , ϕ 1 - k µ ∇P n • ∇C n , ϕ 1 - U n e , ϕ 1 = 0 Cn+1/2 = log(C n+1/2 ) Cn+1 -Cn+1/2 ∆t , ϕ 2 + k µ ∇P n+1 , ∇ϕ 2 + ∂Ω k µ ∂ n P n+1 ϕ 2 = 0 P n+1 -P n ∆t , ϕ 3 + E µ k∇P n+1 , ∇ϕ 3 - ∂Ω E µ k∂ n P n+1 ϕ 3 - EL R e (ψ -P n + RT C n , ϕ 3 - V s U n e , ϕ 3 = 0 C n+1 = exp( Cn+1 ).
This scheme is implemented using the software FreeFem++ [START_REF] Hecht | New Developments in Freefem++[END_REF]. The solution P and C are computed with P 2 and P 1 elements respectively whereas the pure convection terms are solved using the Characteristic-Galerkin method [START_REF] Pironneau | Stability and convergence of a Galerkin-characteristics finite element scheme of lumped mass type[END_REF].

Numerical simulations

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

Validation and comparison with an existing model

With no other 2D dynamic model of phloem transport available, we consider a case analogous to a 1D system so as to compare with the results of [START_REF] Thompson | Application of a Single-solute Non-steady-state Phloem Model to the Study of Long-distance Assimilate Transport[END_REF], a reference implementation of coupled water-solute transport. We simulate phloem transport in a 5 m-long, 10 cm-wide domain with a spatial step ∆x = 5 mm. A 24-hour period is simulated with a time step ∆t = 1 s. The simulation starts with initial pressure and carbohydrate concentration set to zero. Sap viscosity is a function of local solute concentration as a 15 th order polynomial fitted on experimental viscosity for sucrose solutions at T = 293K [START_REF] Swindells | Viscosities of Sucrose Solutions at Various Temperatures: Tables of Recalculated Values[END_REF]. Parameter values are given in Table 1. They are chosen to be equivalent to those used in [START_REF] Thompson | Application of a Single-solute Non-steady-state Phloem Model to the Study of Long-distance Assimilate Transport[END_REF]. Figure 2 shows the predicted profiles of sucrose concentration (C), hydrostatic pressure (P) and axial water flux (J). All profiles are qualitatively very similar to those predicted in [START_REF] Thompson | Application of a Single-solute Non-steady-state Phloem Model to the Study of Long-distance Assimilate Transport[END_REF]. Peak positions, gradients magnitude and transitions over time are well reproduced. The profiles are also quantitatively comparable: the difference is within a few percent. It is only near steady-state, at t = 24 h, that the variation becomes significant for C (less than 10% underestimation) and P (less than 20% underestimation). Although the behaviour at t = 24 h is not particularly physicalsimulated P is twice as high as the highest known measured value (2.4 MPa, [START_REF] Fisher | Sieve tube unloading and post-phloem transport of fluorescent tracers and proteins injected into sieve tubes via severed aphid stylets[END_REF]), loading and leaf dynamics are periodic, not constant [START_REF] Feugier | Dynamical feedback between circadian clock and sucrose availability explains adaptive response of starch metabolism to various photoperiods[END_REF] -both models should yield closer predictions. While the source of the discrepancy has not been identified, it could result from any of the following: differences in parametrisation due to domain geometry, numerical error, missing P derivative in eq. (1.2), implementation details, or the alternate µ = f (C) function.

Symbol

The discrepancy originates from the present model predicting a flow that is slightly more efficient if compared to [START_REF] Thompson | Application of a Single-solute Non-steady-state Phloem Model to the Study of Long-distance Assimilate Transport[END_REF]. As a result less solute accumulates and less pressure builds up in the phloem. The difference is larger near steady-state because error accumulates at each time step. It is also more visible for C and P, which appeared very sensitive to small alterations of the flux. In Figure 3, we show the numerical convergence of the scheme. Here different values of ∆t and ∆x are computed uniformly from 0.5 s to 10 s and from 1 mm to 100 mm respectively. We compute the maximum error by comparison with the approximated solution for ∆t = 0.1 s, ∆x = 0.1 mm after 12 hours. We also plot the lines ∆t and ∆x to demonstrate that the scheme is of order 1.

-1 -0.5

0 0.5 1 1.5 2 2.5 -7 -6 -5 -4 -3 log(∆ t) ∫ Ω C ∫ Ω P O(∆ t) 0 1 2 3 4 5 -7 -6 -5 -4 -3 -2 -1 0 log(∆ x) ∫ Ω C ∫ Ω P O(∆ x)

Parametric study

Several approaches are possible to simplify the governing equations (1.1) and (1.2). Sap viscosity can be treated as a constant instead of a function of sucrose concentration (µ = µ 0 ). This eliminates one term from eq. ( 1.2) and saves the computational cost of evaluating viscosity at all points of the grid. Alternatively, the effect of sucrose's partial molal volume can be neglected (V s = 0). Also, phloem's thickness may also not be updated during the simulations (e = e 0 ). Each approach is considered individually as well as all at once and compared to results of the previous section. Figure 4 shows how the axial water flux is affected by the proposed simplifications. Small changes in the J profile can have large effects on both C and P profiles (not shown here). With a constant sap viscosity, the flux is underestimated during the earlier stages and gets overestimated in the final stage. The position of peak flux also moves more slowly down the conduit with µ = µ 0 . Ignoring the contribution of sucrose to sap volume (V s = 0) causes the flux to be underestimated at any time. Initially, the effect is weak while the sucrose concentration is low but progressively increases in magnitude. Modelling the phloem thickness as being independent of pressure introduces only marginal differences. Although the longitudinal gradient of J is slightly higher between the loading zone and the front of the flow, the magnitude is comparable to that of the reference simulation. Because the simplification has only a minor influence on the J profile, it may seem advantageous to avoid computing phloem's deformation under flow. However, the simplification also has little consequences from a computational point of view; it does not modify the governing equations and all terms must still be evaluated. As expected [START_REF] Lacointe | Modelling phloem and xylem transport within a complex architecture[END_REF], the flux predicted by combining all simplications is very close to the reference behaviour in near steady-state conditions. On the other hand, the flux in the early phases (0:10, 1:10) is strongly affected. Those simplifications are thus best avoided when attempting to describe phloem dynamics and rapid transitioning. Overall, none of the proposed simplifications appeared to be sufficiently beneficial to be included in the model.

Towards realistic designs

The objective of this application is to depict how the model can help in studying the behaviour of living trees. We simulate sap flow on the silhouette of an entire tree. In this application, a finite element mesh is created by Delaunay triangulation from a photograph of Te Matua Ngahere (see fig. 5). This kauri (Agathis australis) tree is the second largest in New Zealand. The height, diameter and volume of the tree trunk are equal to 10.21 m, 5.22 m and 208.1 m 3 , respectively. There are obvious limitations to this approach. The surface is reconstructed from a photographic projection and flat. It means that the mesh is geometrically distorted compared to the actual 3D phloem layer and the results in this application get more inaccurate towards lateral edges. The main objective here is to show that the model can operate on a geometry that is not limited to rectangular channel but that has been retrieved from natural objects. Figure 6 shows the tentative profiles for sucrose concentration and pressure in the phloem. Because of the aforementioned limitations, the profiles are not expected to be realistic. To improve on the quality of the results would require further steps such as reconstructing the surface of an existing tree stem in 3 dimensional space, determining the pattern of sieve cells' orientation on that surface, and identifying the strength of all carbon sources and sinks. Undertaking those steps is beyond the scope of this study.

Orthotropic transport

In this application, the potential for interaction between orthotropic (direction-dependent) transport and competing sinks is investigated. A simulation is carried out for a plate of dimensions L = 1 m, w = 0.6 m and e = 7.5 µm. Three sources, denoted r i ∀i ∈ [START_REF] Thompson | Application of a Single-solute Non-steady-state Phloem Model to the Study of Long-distance Assimilate Transport[END_REF][START_REF] Horwitz | Some simplified mathematical treat-ments of translocation in plants[END_REF], are located near the top of the conduit. They are aligned diagonally to represent the spiral arrangement of branches (phyllotaxis). They have identical strength with a loading flux equal to U . The botton region of the plate (y < 0.1 m) is divided into two sinks, s 1 on the left-hand side (x < 0.3 m) and s 2 on the right-hand side (x > 0.3 m). The unloading at sinks is equal to U i = p i U C C * with p 1 = 0.5 and p 2 = 1. The lateral permeability is set to k T = 0.01k L . The value is set arbitrarily to induce anisotropy in the system but not so high as to completely inhibit lateral transport. It represents the intermediate between isotropic transport (k T = k L ) and the very high ratio (k T = 10 -4 k L ) suggested by [START_REF] Sellier | Phloem transport in trees: A generic surface model[END_REF] to recreate experimental conditions.

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

Three-dimensional surfaces

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

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

Conclusion

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

Major challenges have to be addressed before phloem transport can be computed for real tree configurations. The first challenge is to generate the external surface of an experimental tree, which could be done using 3D laser scanner or terrestrial LiDAR for trees of moderate size. Another challenge is to identify the distribution of diameter and orientation of sieve elements in a tree. The flow grain analogy [START_REF] Phillips | Flow grain analogy[END_REF] or models of auxin transport [START_REF] Kramer | Wood Grain Pattern Formation: A Brief Review[END_REF] could assist in approximating the pattern of fibre orientation. A third, very significant challenge is to evaluate the individual strength of all carbon sinks and sources within the tree over time. The emergence of experimental techniques relying on radiotracers and Positron Emission Tomography [START_REF] Nagel | Temperature responses of roots: impact on growth, root system architecture and implications for phenotyping[END_REF][START_REF] Rotsch | Radiosynthesis of 6-Deoxy-6[18F]Fluorosucrose via Automated Synthesis and Its Utility to Study In Vivo Sucrose Transport in Maize (Zea mays) Leaves[END_REF] will be invaluable to monitor in vivo carbon dynamics and to inform models.

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

1 )Figure 2 :

 12 Figure 2: Sucrose concentration (C), hydrostatic pressure (P) and water flux (J) in the phloem as simulated with the model of [1] (solid line) and with the proposed model (dashed line)

Figure 3 :

 3 Figure 3: Maximum of the relative error obtained by comparison with the approximated solution for ∆t = 1 s, ∆x = 5 mm after 12 hours.

Figure 4 :

 4 Figure4: Effects of model simplications on the axial water flux at an early (10 minutes), intermediate (1h10) and late stage (24h). The considered simplifications are a constant viscosity (µ = µ 0 , dashed line); neglecting the partial molal volume of sucrose (V s = 0, dash-dot line); initial geometry (e = e 0 , dotted line); all simplifications combined (grey, solid line); no simplifications (black, solid line).

Figure 5 :

 5 Figure 5: Mesh reconstruction of Te Matua Ngahere from a photograph (credit: D. Sellier).

Figure 6 :

 6 Figure 6: Phloem pressure (P), sucrose concentration (C) and water flux (J) distributions on a Te Matua Ngahere after 12 hours.

Figure 7 :

 7 Figure 7: Phloem pressure (P), sucrose concentration (C) and water flux (J) distributions on an orthotropic plate after 12 hours.

Figure 8 :

 8 Figure 8: Phloem pressure (P), sucrose concentration (C) and water flux (J) distributions in a 3D fork after 12 hours.

Table 1 :

 1 Description of parameters employed in the model. Numerical values correspond to the initial values used in section 3.1.

		Description	Value	Units
	L R	Radial hydraulic conductivity	1.57 × 10 -13	m Pa -1 s -1
	ψ	Xylem hydrostatic pressure	0	Pa
	R	Gas constant	8.31	J mol -1 K -1
	T	Temperature	293	K
	V s	Partial molal volume of sucrose	2.15 × 10 -4	m 3 mol -1
	k L	Longitudinal permeability	9.28 × 10 -12	m 2
	k T	Tangential permeability	9.28 × 10 -13	m 2
	e	Phloem thickness	7.5 × 10 -6	m
	E	Phloem Young's modulus	1.7 × 10 7	Pa
	U	Loading rate	3.375 × 10 -6	mol m -2 s -1
	C *	Reference sucrose concentration	500	mol m -3
	µ	Viscosity	10 -3 at C = 0	Pa s
	L R RT -Vs U C *	Positivity constraint	3.808 × 10 -10

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