Jean Claude Mailhol

Pierre Ruelle

Sabine Walser

Niels Schütze

Cyril Dejean

Analysis of AET and yield predictions under surface and buried drip irrigation systems using the Crop Model PILOTE and Hydrus-2D

Keywords: Actual evapotranspiration, crop model, drip irrigation, Hydrus-2D

Innovative irrigation solutions have to face water scarcity problems affecting the Mediterranean countries. Generally, surface (DI) or subsurface drip irrigation systems (SDI) have the ability to increase water productivity (WP). But the question about their possible utilization for crops such as corn would merit to be analysed using an appropriate economic tool. The latter would be necessary based on the utilisation of a modelling approach to identify the optimal irrigation strategy associating a water amount with a crop yield (Yc). In this perspective, a possible utilisation of the operative 1D crop model PILOTE for simulating actual evapotranspiration (AET) and yield under a 2D soil water transfer process characterizing DI and SDI was analysed. In this study, limited to a loamy soil cultivated with corn, the pertinence of the root water uptake model used in the numerical code Hydrus-2D for AET estimations of actual evapotranspiration (AET) under water stress conditions is discussed throughout the Yc=F(AET) relationship established by PILOTE on the basis of validated simulations. The conclusions of this work are (i): with slight adaptations, PILOTE can provide reliable WP estimations associated to irrigation strategies under DI and SDI, (ii): the current Hydrus-2D version used in this study underestimates AET, compared with PILOTE, in a range varying from 7% under moderate water stress conditions to 14% under severe ones,(iii): A lateral spacing of 1.6 m for the irrigation of corn with a SDI system is an appropriate solution on a loamy soil under a Mediterranean climate.

A local Yc= F(AET) relationship associated with a Hydrus-2D version taking into account the compensating root uptake process could result in an interesting tool to help identify the optimal irrigation system design under different soil conditions.

Introduction

Mediterranean countries are more and more confronted to water scarcity problems. Therefore, farmers are encouraged to adopt efficient irrigation systems. That is the case yet in most part of the Maghreb countries with the intensification of drip irrigation. When the water resource must be shared with other users, supplementary efforts from farmers can be required to save water by optimizing their irrigation strategies. Modelling could be the adequate tool for conducting this optimization where a crop model and an irrigation technique would be narrowly linked. Consequently, analysing the impact of an irrigation strategy on the crop yield, whatever the irrigation system, appears as a necessity. Although high frequency is generally touted as a major advantage of micro irrigation, studies have shown that it is not generalizable depending on soil, crop type and its environmental context as for corn for instance [START_REF] Lam | Subsurface drip irrigation for corn production: a review of 10 years of research in Kansas[END_REF].

Numerical simulations are efficient tools to investigate optimal drip management practices and system design [START_REF] Meshkat | Modeling of evaporation reduction in drip irrigation system[END_REF][START_REF] Assouline | The Effects of Microdrip and Conventional Drip Irrigation on Water Distribution and Uptake[END_REF][START_REF] Schmitz | New strategy for optimizing water application under trickle irrigation[END_REF][START_REF] Cote | Analysis of soil wetting and solute transport in sub-surface trickle irrigation[END_REF][START_REF] Skaggs | Comparison of Hydrus-2D simulations of drip irrigation with experimental observations[END_REF][START_REF] Beggs | Modeling subsurface drip application of onsite wastewater treatment system effluent[END_REF][START_REF] Li | Modeling of water flow and nitrate transport under surface drip fertigation[END_REF][START_REF] Lazarovitch | System dependent boundary condition for water flow from subsurface source[END_REF]2007;[START_REF] Gardenas | Two-dimensional modeling of nitrate leaching for various fertigation scenarios under micro-irrigation[END_REF][START_REF] Hanson | Evaluation of urea-ammonium-nitrate fertigation with drip irrigation using numerical modeling[END_REF]2008). But the latter refer only to the water transfer aspect and do not account for the impact of an irrigation strategy on the water plant uptake and on the more reason on crop yield. Coupling a soil water transfer model with a crop model would serve adopting the best design system with its irrigation strategy adapted to the water constraints on the basis of economical consideration (farmer acceptability).

Few researches refer to that coupling problem. It was addressed by Mmolawa and Or (2003) when testing the analytical model of [START_REF] Coelho | A parametric model for two-dimensional water uptake intensity[END_REF] for soil water dynamic with spatially plant water uptake based on a parametric model for the root system pattern and the model of [START_REF] Warrick | Time dependent linearised infiltration . I. Point sources[END_REF] for the water transfer. The latter which assumes that the derivative of hydraulic conductivity by soil water content is constant, is often criticized. This assumption is considered as very strong by [START_REF] Ababou | Modélisation des transferts hydriques dans le sol en irrigation localisée PhD Univ[END_REF], and gives results too far from observations according to [START_REF] Clothier | Constant -flux infiltration from a hemispherical cavity[END_REF]. In their article Mmolawa and Or (2003) recognized improvements which are still required for a better prediction of the water dynamic within the simulated domain during irrigation but also to better account for the plant water uptake process during redistribution. Coupling a crop model with a 2D soil water transfer model such as Hydrus-2D [START_REF] Simunek | The Hydrus-2D software package for simulating the two-dimensional movement of water, heat, and multiple solutes in variably-saturated media[END_REF] could be realized for the case of a linear source (at the soil surface or buried) as done for the furrow case by Wöhling et al., (2007) and [START_REF] Wöhling | A physically based model for simulating 1D surface-2D subsurface flow and plant water uptake in irrigation furrows. II Model test and evaluation[END_REF]. But it is of prior interest to analyse the reliability of some modelling alternatives much more operative through a rigorous experimental protocol to highlight the complexity of processes involved in the Soil-Plant -Atmosphere system under drip irrigation. This analysis is encouraged by the following developments which still emphasize the choice of a lumped model for dealing with the plant water uptake phenomenon when focusing on AET and yield prediction only. Indeed, basing plant water uptake from a map of the root density only, would consist to ignore the capacity of the plant to remove its water uptake capacities to soil regions where water is more available when a restrictive irrigation strategy is applied. The case of drip irrigation even more emphasizes the root water uptake problem under water stress conditions because plants not only remove their water uptake capabilities but can also remove sensitively their root system in the vicinity of the source (Mubarak et al., 2009b).

For unstressed conditions root water uptake models based on a observed root density pattern as proposed in Mubarak et al. (2009b) or a parametric one (Coehelo and Or (1996), are less or more efficient, the reference in term of water uptake being still that implemented in Hydrus-2D However, these models do not integrate yet the fact that under water stress conditions most of the water is absorbed from the zone of low tensions. Thus a relative small part of the root system can be responsible for most of the plant water uptake [START_REF] Feddes | Séminaires sur l'irrigation localisée[END_REF]. Water is taken up by roots from various soil depths in a manner which maximizes the plant's chances for survival [START_REF] Lubana | Modeling soil water dynamics under trickle emitters-a review[END_REF]. From their sound literature analysis of plant water uptake models, [START_REF] Yadav | Soil moisture dynamics modelling considering the root compensation mechanism for water uptake by plants[END_REF] showed the difficulty to satisfactory model such complex process mixing both physical and physiological factors. The results of their research show that under favourable soil moisture conditions, plants extract water at the maximum rate according to the root distribution pattern and when the moisture stress occurs in the upper soil profile, the diminished water uptake from the lower water scarce region is compensated by an enhanced water uptake from lower wetter layers. Aware of this problem, [START_REF] Simunek | Compensated Root Water and Nutrient Uptake Ecological modelling[END_REF] proposed a calculation procedure to mimicking this physiological behaviour allowing the plant to increase its transpiration rate. This new formulation of the plant water uptake could be taken into account in a new version of the Hydrus-1D and 2D codes. But the efficiency of this new formulation in the AET estimation by Hydrus-2D will require that root growth can be simulated, which is not the case yet.

The objective of this paper is to show that a crop model such as PILOTE [START_REF] Mailhol | AET and yields assessments based on the LAI simulation. Application to sorghum and sunflower crops[END_REF][START_REF] Khaledian | Adapting PILOTE model for water and yield management under direct seeding system (DSM). The case of corn and durum wheat in mediterranean climate[END_REF], a crop model presumed limited to simulations under sprinkler or under border irrigation [START_REF] Mailhol | SPFC : a tool to improve water management and hay production in the Crau region[END_REF], can be used for AET and yield predictions for a surface and buried drip irrigation system meaning some utilisation conditions or model adaptations A sound analysis of the water transfer under DI and SDI will be carried out under cover of a rigorous experimental protocol for attempt understand the fate of water under such irrigation systems applied to field crops. This analysis will be completed by a modelling approach. The latter will be based on the conjoined utilisation of PILOTE and Hydrus-2D, a model theoretically more adapted to a spatially variable soil wetting pattern. Field results and modelling analysis should contribute to a better understanding of the plant water uptake process under restrictive irrigation conditions.

Furthermore, the paper proves that SDI can improve water productivity of corn under a mediterranean climate.

Materials and methods

The experimental setup and measurements

The experiments related to three irrigation seasons (2007,2008,2009). They were conducted at the experimental station of Lavalette at the Cemagref Institute, Montpellier, France (43°40'N, 3°50'E) on a loamy soil plot with corn. The soil properties of this loamy soil are presented in Table 1. A fully equipped meteorological station exits on the site. Maize (Pionner PR35Y65) was sown on April 24 for the three seasons excepted for a specific treatment of 2009. The experimental protocol of 2007 is described in Mubarak et al., (2009). It is recalled here with regards to the specificity of the study. Corn was sown with a row spacing of 0.75 m and with a density of 100000 plts ha -1 . The surface drip irrigation system (DI) was installed at the beginning of June with one drip tube for two plant rows resulting in lateral spacing of 1.5 m. The tubing was a drip tape with a 16-mm inside diameter, a wall thickness of 8 mm, an emitter spacing of 30 cm and a flow rate of 3.67 L h -1 m -1 at a pressure of 0.8 bars. First irrigation was applied at the 12-leaf stage of development of the crop about two months after sowing i.e., on June 19, for a period of approximately 7.5 h. Due to rainy events it was not necessary to irrigate early as attested in Mubarak et al. (2009a. On the basis of estimated crop water demand, irrigation was applied daily for three hours in a full irrigated treatment (FIT) and every other day in a limited treatment (LT). As it was not possible to irrigate at the weekend, three hours supplementary irrigation was applied on Friday and Monday. This irrigation schedule was maintained until the last week of July. Then, after an interruption of few days, irrigation started again but at a rate sensitively modified compared with the previous one.

For each treatment, two access tubes for the neutron probes were installed just before irrigation started to a depth of 150 cm, the measurements being done at a step depth of 10 cm. One was placed beside of the central drip line and the other in the crop row, with 37.5 cm between the two. The neutron probe was calibrated from gravimetric measurements using soil samples collected during installation of the access tubes.

Leaf area index (LAI) was monitored along the cropping cycle, measurements being done approximately each week using the LAI2000 (LI-COR Plant Canopy Analyzer). Plant samples are collected after maturity for evaluating dry matter and grain yield according to the protocol described in [START_REF] Mailhol | AET and yields assessments based on the LAI simulation. Application to sorghum and sunflower crops[END_REF] and [START_REF] Khaledian | Adapting PILOTE model for water and yield management under direct seeding system (DSM). The case of corn and durum wheat in mediterranean climate[END_REF]. At least, just after harvesting, a root density profile has been established as described by Mubarak et al., (2009b).

In 2008, the DI system was completed by a subsurface drip irrigation system (SDI) installed at a 35 cm depth. It was used with two different lateral spacing of 160 and 120 cm respectively and centered between the corn rows spaced from 80 and 60 cm respectively. The whole resulted in two irrigated treatments, SDI160 and SDI120 of 1200 m 2 (60 m length, 20 m width).Although having a holding capacity close to 180 mm/m, the clay content of these two treatments is slightly different, resulting in a volumetric soil water content at field capacity of 0.3 for SDI160 and 0.32 for SDI120, soil water content at wilting point being taken at 0.12 and 0.14 respectively. The buried pipes were interconnected at their ends by a stabilization pipe. The latter was equipped with a regulation system, a drain-cock and manometers to check the conditions of functioning. Volumetric valves allowed a precise knowledge of the water application depths The flow rate was q = 2.8 L h -1 m -1 in average at a pressure of 0.7 bars, values lower than for DI. Plant density was 100000 plts.ha -1 in SDI120 due to a row number higher than in SDI160 where plant density was 78000 plts.ha -1 only. Irrigation was applied three times a week and the dose was doubled on Monday and Friday, due to the fact that there was no irrigation during the week-end. It was established in order to approximately meet 70% of the maximum water requirements, the objective being also to deliver the same water amount for the two treatments. The total water amount delivered to the two treatments was of 230 mm.

For each treatment an access tube for the neutron probes was installed at the vertical of the crop row. For analyzing the evolution of the wetting pattern under irrigation, two capacitance sensors place in access tube (EnviroSMART) were vertically installed, a first under row and, a second at 15 cm from the buried dripper line. Note that fertigation was used for surface and subsurface drip irrigation. Nitrogen applications were adjusted to satisfy the plant requirements on the basis of the initial soil nitrogen content [START_REF] Khaledian | Adapting PILOTE model for water and yield management under direct seeding system (DSM). The case of corn and durum wheat in mediterranean climate[END_REF].

At harvest, a soil profile was dug by a backhoe to analyze the rooting system of the crop perpendicularly to the drip lines for SDI160 and SDI120. Root repartition was described on a grid of 5x5 cm (Fig. 1 ) using a qualitative root density index, described by Mubarak et al., (2009b) that varies from 0 to 5 when respectively there no roots or the root quantity is maximum.

In 2009, two rainfed treatments (zero irrigation) RF60 (row spacing = 60 cm) sown on 05/07 and RF75 (row spacing = 75 cm) sown on 04/23, were monitored. A replication of the 2008 subsurface system was implemented but with SDI120 sown 15 days latter (on 05/07) than SDI160. But unfortunately, due to problems during seeding, the standard configuration of 2008 assuming that a lateral is centered between two crop lines is far to be respected. Indeed, the row spacing is of 0.75 instead of 0.8 on SDI160 and of 0.60 for SDI120. The resulted arrangement of the crop rows and drip lines was non-uniform and consequently, some crop rows were remote to irrigation drip lines when others were very close.

The soil water content monitoring in 2009 was improved by the installation of supplementary access tubes for the neutron probes. In addition to that installed vertically to the crop row, a second one was installed near a buried tube (Fig. 2a), and, only for SDI160, a third one between two rows (Fig. 2b). As shown in Fig. 2b for SDI160, the distance between the access tube under crop and the lateral is not so far as presumed. For the two SDI treatments, a series of mercury tensiometers were installed at 10,20,30,45,60,75,90,110,130, and 150 cm depth on a crop line in the vicinity of the neutron probe access tube. As in 2008 two capacitance sensors (EnviroSMART) were vertically installed as presented in Fig. 2d.

The flow rate was approximately the same as in 2008, but the irrigation scheduling of 2009 strongly differed from that of 2008. The growth period can be characterized as dry. As shown on Fig. 3a, rain mainly stopped at the end of April in 2009, interrupted by little rainfall in May (18 mm on 07/05) and in June (10 mm on 06/07) while it continued until mid June in 2007 and 2008. A different irrigation scheduling was tested in 2009. It consisted of irrigating at low frequency but with high water application depths, the irrigation duration overreaching sometimes 10 hours resulting in water application depths (WAD) greater that 30 mm as presented in Fig. 3b.This irrigation schedule was obtained by a stochastic SDI management framework [START_REF] Walser | Towards higher crop water productivity using a simulation based controlled deficit irrigation strategy[END_REF]. The framework developed to realize a real-time scheduling is based on a stochastic approach which uses a set of climate series consisting of observed weather data, actual weather forecasts and weather scenarios which are representative for the long term climate pattern until the end of the growing season. The simulation/optimization package comprising PILOTE and GET-OPTIS [START_REF] Schütze | Novel simulation-based algorithms for optimal open-loop and closed-loop scheduling of deficit irrigation systems[END_REF][START_REF] Schmitz | GAIN-P A new strategy to increase furrow irrigation efficiency[END_REF] generates then a set of corresponding schedules which are evaluated by a statistical analysis in order to provide a reliable schedule for the actual date of the growing season. The LARS-WG weather generator [START_REF] Semenov | Comparison of the WGEN and LARS-WG stochastic weather generators for diverse climates[END_REF] was used to analyse the long term climate characteristics based on Cemagref historical climate data from 1991-2008 and to provide representative weather scenarios. The SDI management framework was run once a week in order to adapt the irrigation schedule to the actual weather.A water application depth (WAD) in mm is obtained by multiplying q(L h -1 m) by the irrigation duration and divided by the pipe spacing. The irrigation dose was distributed on two days when the tensiometer readings indicated drainage risks.

The protocol of the plant sample collection for yield measurements has been modified to account for the visible yield heterogeneity on SDI160 and SDI120 resulting from the fact that a plant row can be near from the pipe (the N row) and the adjacent one can be far (the S row). Consequently a yield sub-plot was established by the two rows N and S. On SDI160 the plants were harvested on a length varying from 2m (be 10 plants in average collected by row) the sub-plot width being of 1.5 m. On SDI120 plants were harvested on a row length of 2.5 m (be 15 plants in average by row), the sub-plot width being of 1.3 m. Plants collected on N et S were weighted separately for constituting a sub-sample assigned to an areas S/2 for a given sub-plot of area S.

Modelling

Water transfer under drip irrigation is axially symmetrical when considering an individual dripper. For a linear source, with a distance between two emitters is of 30 cm water transfer can be assumed to be ruled by a 2D process [START_REF] Skaggs | Comparison of Hydrus-2D simulations of drip irrigation with experimental observations[END_REF][START_REF] Patel | Dynamics and modeling of soil water under subsurface drip irrigated onion[END_REF]Mubarak et al., 2009b). This multidirectional water transfer rise questions about the interest of a 1D model such as PILOTE for DI or SDI.

As previously evoked, PILOTE applications to SDI are here limited to AET and yield simulations since a narrow relationship exits between these two factors [START_REF] Howell | Relationship of dry matter production of field crops to water consumption[END_REF][START_REF] Jones | Plant and microclimate[END_REF] what is confirmed by which of Fig. 4 derived from PILOTE simulations for the two corn varieties studied at Lavalette [START_REF] Khaledian | Adapting PILOTE model for water and yield management under direct seeding system (DSM). The case of corn and durum wheat in mediterranean climate[END_REF] between 1998 and 2007 irrigated with different irrigation systems (sprinkler and DI). LAI and the soil water reserve (SWR) evolution are amongst the main PILOTE outputs. SWR on the maximal root depth Px, will be widely analysed on the basis of a comparison with those calculated from the soil water content profiles measured under crop and under pipe.

For the specific case of SDI, PILOTE needs to be adapted, while for DI, the reduction of soil evaporation (only a part of the soil surface being wetted) can be performed by the model option proposed in [START_REF] Khaledian | Adapting PILOTE model for water and yield management under direct seeding system (DSM). The case of corn and durum wheat in mediterranean climate[END_REF] if necessary (when irrigation starts with LAI <<3).

Hydrus-2D is particularly well adapted to this 2D water transfer. It is used in this study mainly for dealing with the problem of the plant water uptake process under severe water stress conditions. In Hydrus-2D, SWR will be simulated both under crop and under pipe at the opposite of PILOTE which can only give a lumped formulation of the soil water balance illustrated by the evolution of SWR. At last, AET(Hydrus) calculated on the basis of a spatially distributed wetting pattern will be compared to AET(PILOTE).

2.2.1. Adapting PILOTE to SDI PILOTE simulates soil water balance and crop yield at a daily time step by association of a soil module and a crop module, under the assumption of water being the only limited condition. The soil module consists of a three-reservoir system [START_REF] Mailhol | Operative modeling for evapotranspiration assessment : Calibration methodology[END_REF][START_REF] Mailhol | AET and yields assessments based on the LAI simulation. Application to sorghum and sunflower crops[END_REF]) covering surface layer until the maximum rooting depth. The reservoir with shallow depth of 10 cm rules the water balance at the soil surface, in which evaporation is governed by current LAI acting on the partitioning coefficient between transpiration and evaporation. The following reservoir R 2 accounts for root section, so its capacity increases with root growth. Before the potential root area is totally taken by the second reservoir, the third reservoir represents the remaining part. Water is first taken from the shallow reservoir until total depletion by evaporation and plant then, from the second one by plant only. On the basis of field capacity (Fc) and wilting point (Wp), the soil water balance among reservoirs is thus calculated. Maximal evapotranspiratin (MET), and actual evapotranpiratoin (AET) are involved in the water stress index (WSI) calculation. MET is derived from MET = Kc.ETref, where ETref is the reference evapotranspiration [START_REF] Allen | An update for the calculation of the reference evapotranspiration[END_REF] and Kc, the crop coefficient as a function of LAI. Under water stress conditions, AET linearly decreases from MET with the depletion level of R 2 . Then, WSI, obtained accordingly to this lumped plant uptake approach, is exported to the crop module as environment coefficient.

The crop module is based on the LAI simulation and its response to WSI. The simulation involves two shape parameters and a vegetative stage parameter (T m ) corresponding to the temperature sum when the maximum LAI (LAI X ) reached. T m and LAI X can be derived from the literature or measured in the field. Total dry mater (TDM) is calculated based on Beer's Law, RUE (the radiation use efficiency) being affected by WSI. Grain yield is evaluated by the product of TDM by a harvest index (HI). HI is set to a potential value if average LAI from the stage "grain filling" to the stage of '"pasty grain" is greater than a threshold value, otherwise it linearly decreases [START_REF] Mailhol | Analysis of irrigation systems and irrigation strategies for durum wheat in Tunisia[END_REF][START_REF] Khaledian | Adapting PILOTE model for water and yield management under direct seeding system (DSM). The case of corn and durum wheat in mediterranean climate[END_REF]. The required climatic data are precipitations, global radiation, average temperature and ETref..

In this adaptation to SDI, the irrigation dose of the day j: Dose(j) is added to d' 1 (j) the water that drains from R 1 when the later overreaches its field capacity:

d 1 (j) = d' 1 (j) + Dose (j)
(1),

The eventual percolation water d' 1 (j) from R 1 is thus completed by Dose (j), resulting in d 1 (j) which supplies R 2 . This water supplied to R 2 is thus not affected by evaporation but transpired only or eventually partially drained.

The utilisation of such a model under DI or SDI conditions assumes that all the water applied is stored within the root reservoir. This can be realized when no macro pore effect, initiating preferential pass flows, exist or when the irrigation strategy minimizes the drainage risks. Note that due to the spatially variable soil wetting pattern under drip irrigation, drainage can locally occur under the drip line when irrigation is frequently applied even at a rate adequate to plant transpiration [START_REF] Hanson | Leaching with subsurface drip irrigation under saline, shallow groundwater conditions[END_REF].

Modelling with Hydrus-2D

The ability of Hydrus-2D for simulating the water transfer process under drip irrigation has been widely demonstrated on bare and cultivated soils [START_REF] Vrugt | Calibration of a two dimensional root water uptake model[END_REF][START_REF] Skaggs | Comparison of Hydrus-2D simulations of drip irrigation with experimental observations[END_REF][START_REF] Assouline | Soil-Plant System Response to Pulsed Drip Irrigation and Salinity: Bell Pepper. Case Study[END_REF][START_REF] Arbat | Monitoring soil water status for micro-irrigation management versus modelling approach[END_REF][START_REF] Patel | Dynamics and modeling of soil water under subsurface drip irrigated onion[END_REF]. The ability of this model for simulating the water transfer process for the corn case at Lavalette has been widely demonstrated in [START_REF] Mailhol | Impact of fertilisation practices on nitrogen leaching under irrigation[END_REF], [START_REF] Wöhling | A physically based model for simulating 1D surface-2D subsurface flow and plant water uptake in irrigation furrows. II Model test and evaluation[END_REF] for furrow irrigation, and by Mubarak et al., (2009b) for DI. The application of this model is only used here for the SDI system for evaluating the main component of the soil water balance under this SDI system with regards to the irrigation strategy (high WAD at low frequency), and the evoked asymmetric problem. This simulation results will be then compared with those evaluated by PILOTE.

The simulated domain is comparable to which presented in Mubarak et al., (2009b). The width of the domain corresponds to the distance separating the two zero flux plans: 1.6 m for SDI160 for instance. But here, the drip tubing is buried at 35 cm depth. The configuration of the root pattern adopted in Hydrus-2D differs from that adopted by Mubarak et al. (2009b) to account for the fact that one of the two plants was close to the pipe when the other was remote. The constitution of that uptake pattern is based on the root density measurements realized in the vicinity of the neutron probe sites by the same method as which described for 2008.

The soil hydraulic parameters are drawn from Mubarak et al., (2009b) and presented in Table 1. The same procedure as [START_REF] Skaggs | Comparison of Hydrus-2D simulations of drip irrigation with experimental observations[END_REF], [START_REF] Patel | Dynamics and modeling of soil water under subsurface drip irrigated onion[END_REF] and Mubarak et al., (2009b) was adopted for calculating the drip line boundary water flux q (m h -1 ). The latter is multiplied by the daily water application duration because our objective was mainly to simulate the main water balance components along the cropping cycle with a reasonable calculation time. Thus, the atmospheric conditions were taken into account on a daily basis. Crop water demand ETc was calculated by multiplying the ET ref , by the crop coefficient Kc linked to LAI, as in PILOTE, with a maximal value of Kc (Kc max = 1.2) proposed by Doorenbos and Pruit, (1977). As Hydrus-2D does not simulate root growth, the simulation period is divided into two periods, the soil water content profile at the end of the first period being used for the initial soil water content profile of the second period. The root density for each period is obtained as proposed in Mubarak et al., (2009b). Nevertheless, one must point out the difficulty of proposing an exact root pattern evolution for different dates along the cropping cycle on the basis of that established just after harvesting. Depending on the adopted root pattern for a given period, the soil water content pattern within the simulated domain can be significantly different. A solution for proposing an adequate root pattern corresponding to a simulated period can be obtained by correcting the Hydrus-2D root pattern from neutron probe profiles using an errors and trial approach. Under no water stress conditions, that manner of proceeding can help identify the real root pattern. But, under water stress conditions, the identified root pattern could be that obtained to mimic the compensation root water uptake process. For instance, by transferring the plant water uptake capabilities of the upper layers to the deeper layers it is possible to improve the soil water content simulations of the deeper layers for water stress periods.

Experimental results

3.1 The root system under SDI A mean value of the root index was calculated for each depth and on the North side and South side of the profile to examine possible effect of drip tape and its wet bulb on root repartition. The value of that root index as a function of the depth is presented in Fig. 5 for the SDI160 in 2008 (a) and in 2009 (b) in the asymmetric case (pipe almost under a crop row). Maximum rooting depth is around 1.2 m and most of the roots are observed in the first 40 to 50 cm for all cases. In the asymmetric case, root density is much higher under pipe (35 cm at the vertical of the row) than in 2008 where the pipe is between the rows.

In Fig. 6 the map of the root index (RI) indicates that the roots are present in the whole of the domain even far from the buried pipe. High RI values often coincide with the local presence organic matter (crop residues). That still underlines the existing difficulty to model the root system distribution for simulating the plant water uptake process. For instance, the presence of roots more especially in the vicinity of the soil surface does not mean that they remained active along the cropping cycle as further seen. The high RI values close to soil surface is probably due to the fact that rainy events at the end of April (33 mm) and at mid May (17 mm) have contributed to a good installation of the root system of the SDI160 treatment. A same configuration of RI can be noticed for SDI120, but with lower RI values, maximal RI not exciding 3.5, due to a later sowing date than SDI160. After observations done on other sites, when the drip line is under a corn row, (5 cm apart as shown on Fig. 6,) a high increase of root density can be seen under the drip line for a depth of 40 cm or more. A same effect existed for SD120 with a drip line 15 cm apart the corn row. It can be concluded that drip line has a noticeable effect on root development when the distance drip line-crop is short. Such root density configuration was not observed in 2008 with laterals placed at mid distance from crop rows.

Analysis of soil water content profile from neutron probe measurements

In 2008, soil contribution to plant water supply is particularly well highlighted on SDI160 confirming that soil participation to the plant water supply legitimates strategies at low irrigation frequencies such as those applied in 2009. For both treatments corn took up water until 1.20 at least while tensiometer monitoring for both SDI treatment in 2008 do not reveal drainage losses.

In 2009, water storages calculated on a thickness of 30 cm, from the soil surface to 150 cm, were continuously decreasing on T1 and T2 for both treatments and are similar under 120 cm at the end of the cropping cycle. The fact that a comparable trend is not observable for T3 of SDI160 results from the lateral spacing, a spacing of 1.2 m only, would allow the maintain of a root activity much greater in the vicinity of SDI120 T2 than in the vicinity of SDI160 T3. On SDI120, under deficit irrigation, crop is also more able to use a high amount of the soil water than on SDI160.

Tensiometer monitoring for both SDI treatments in 2009 revealed drainage losses under pipe (case where the pipe is under crop) during short periods corresponding to a WAD > 20 mm delivered at the beginning of the irrigation season when the plant water uptake rate is low yet.

Analysis of the Soil moisture evolution from EnviroSMART sensors SDI160 case in 2008

From the analysis of the sensor responses to irrigation events one can say that irrigation water reached only a depth of 45 cm: and there is no increase of the signal of the deeper layers (55 and 65 cm). The highest sensor response was obtained at 33 cm depth. For the deeper layers, there was a small decrease of the signal meaning that crop took up water from soil storage. From the analysis of sensor responses installed under the corn row one can say that lateral diffusion did not reach the corn row, meaning that lateral diffusion is lower than 30 cm. SDI120 in 2009 (from June 22 to July 02): an asymmetric case (row at 45 cm from pipe instead of 30 cm for a symmetric case) The analysis focused on SDI120 only. The irrigation effect is not perceptible at the vertical of the crop between 0.1 to 0.6 m depth. That seems logical with such high distances from the dripper line: 45 cm for this asymmetric configuration. Between rows (13 cm from pipe), the water front reached 1 m depth but does not seem to go further. The first irrigation (11.6 mm) did not inverse the decreasing trend of the soil water content under row at the opposite of the inter-row where soil water content variations can be seen down to 0,5 m. The second Author-produced version of the article published in Agricultural Water Management, 2011, 98(6), 1033-1044. The original publication is available at http://www.elsevier.com/locate/agwat irrigation stops the decreasing trend under row and even generated a sensitive increase of the soil water content at a depth of 0.7 m at the opposite of 2008 due to high WADs delivered at a low frequency. That is interesting in term of horizontal water diffusion even at a depth greater than 0.65 m, when assuming this water amount will be profitable to plant and will not contribute to eventual drainage losses. Indeed, the water diffusion process, enhanced by the soil water deficit before irrigation [START_REF] Philip | Travel times from buried and surface infiltration point source[END_REF], would reach parts of the root system at deeper layers. Note that the Ks decrease (Table 1) between first (0-55 cm) and 2 nd layer (55-90 cm) also contributed to this diffusion process as further confirmed by Hydrus-2D. Both effects resulted in a wetting zone display at around 0.7 m depth, the lateral extension of the wetted zone being maximum around 45 cm from the drip line at 70 cm depth

Yield results in 2009

Due to asymmetric problems, the yields present a high variability as attested by Fig. 7, the latter being a little higher for SDI160 than for SDI120. The averaged deviation between yield row N and yield row S is the lowest (3.6 vs 5 Mg/ha) for this treatment. For SDI120, averaged yields are 19 Mg/ha and 11.8 Mg/ha for total dry mater (TDM) and grain yield (GY at 15% humidity) with a coefficient of variation Cv =21%. These values are of 20 and 11.8 Mg/ha respectively for TDM and GY for a same water amount of 250 mm. Note that averaged yield for SDI120 would have probably been higher than the one of SDI160 if both had been sown at the same date. As confirmed by the two rainfed treatments the yields are higher for corn sown earlier: 8.4 (for RF75) vs 7.5 T/ha (for RF60) for TDM and 3.3 Vs 2.5 T/ha for GY. For the full irrigated treatment (the surface drip system) TDM and GY are 25 and 16 Mg/ha respectively for a water amount of 350 mm. For these three treatments, the Cv value is lower than 10%.

PILOTE simulations

The ability of PILOTE for simulating the soil water reserve (SWR) evolution on the maximal root depth Px, LAI and the yields for the corn case was already proved in [START_REF] Khaledian | Adapting PILOTE model for water and yield management under direct seeding system (DSM). The case of corn and durum wheat in mediterranean climate[END_REF] for all the treatments of 2007 sprinkler irrigated. In the frame of this work only the DI results are presented for 2007. As shown by Fig. 8(a) it is remarkable that SWR of FIT simulated by PILOTE, matches well with the mean water storage on Px =1.2 m calculated by ,

  Px dz z SWR 0 ) (



, from the neutron probe access tube installed under crop and which installed at the vertical of the surface pipe A similar statement can be done from Fig. 8(b) regarding the LT treatment where only data of the tube installed under pipe was available until the end of the irrigation period. The decreasing trend simulated by PILOTE, here more perceptible than on FIT highlights the ability of the model for simulating the lumped plant uptake process.

As irrigation started when LAI was greater than 2, the model option allowing a soil evaporation reduction was not activated. The yields are well simulated by PILOTE for these two treatments: 29.2 vs 29.4 and 17.2 vs 17.4 for DM (Mg/ha) and grain yield (at a humidity of 15%) respectively for FIT and 28 vs 27.5 and 16.5 vs 16.4 for LT, the water application depth (WAD) of these two treatments being of 430 and 306 mm respectively. Note that HIpot, the potential harvest index, is set at 0.5 for a plant density of 10 plts/m 2 as in [START_REF] Khaledian | Adapting PILOTE model for water and yield management under direct seeding system (DSM). The case of corn and durum wheat in mediterranean climate[END_REF].

For the SDI treatments of 2008, SWR(PILOTE) does not differ a lot from which measured under crop, when it was expected a much more gaps since PILOTE is based on a 1D water transfer process. Meanwhile, the yields are correctly predicted: 24 vs 24 and 14.7 vs 15 for DM (Mg/ha) and grain yield respectively for SDI160 and 25.6 vs 26 and 15.2 vs 15 for SDI120.

SWR simulations of 2009 on the SDI treatments are presented at Fig. 9. At the beginning of irrigation (from DOY 165 to DOY 180), SWR (PILOTE) under crop and SWR measured under pipe differ a lot at this period. This is due to high and frequent WADs, crop demand being low yet and the soil water reserve still far from depletion. Then they tend to converge because the pipe is close to plant row and the water applications are less frequent and lower. The SWRs are almost confounded, more especially for SDI120 (Fig 9b) at the end of the cropping cycle due to the soil water redistribution process fostered by irrigations at low frequency, SWR(PILOTE) being logically sensitively above measured SWR's. Regarding Fig. 9a, with pipes assumed to be exactly between the corn rows, one would expect to see SWR(PILOTE) close to the averaged SWR's as for DI in 2007. However, with DI, the 2D water transfer process begins from the sol surface. That can explain why SWR(PILOTE) is exactly between SWR(measured under pipe) and SWR(measured under crop). This is not the case with SDI because the water transfer process begins only from about 40 cm, a depth representing 1/3 of the maximal depth on which SWR is calculated. As we shall see further with the case of SDI in 2008, the 2009 asymmetric problem is not the main reason why SWR(PILOTE) does not fall exactly between SWR(under Crop) and SWR(under Pipe). It can be noticed that from DOY210, SWR(inter C-L) remains quite constant, plant probably enhancing its water uptake capabilities where water is much more available i.e. at deeper layers under drippers. That is in agreement with the results of [START_REF] Yadav | Soil moisture dynamics modelling considering the root compensation mechanism for water uptake by plants[END_REF] illustrating the dynamic root compensation mechanism. We nevertheless have to note that the activity radius of the neutron probe (about 15 cm) tends also to mitigate the SWR differences. As for DI, SWR(PILOTE) under SDI follows well enough the measured decreasing trend of SWR, meaning that a lumped water uptake is simulated by PILOTE.

For obtaining reliable AET simulation results, it is necessary that the LAI is simulated correctly since in PILOTE, Kc is derived from LAI., PILOTE simulates well enough the LAI evolution, the coefficient of efficiency (Ce) of Nash-Sutcliffe (ASCE, 1993) being of 0.989 and of 0.941 for SDI160 and SDI120 respectively. The simulation results show that PILOTE slightly better simulates LAI under a sprinkler system than under SDI (see [START_REF] Khaledian | Adapting PILOTE model for water and yield management under direct seeding system (DSM). The case of corn and durum wheat in mediterranean climate[END_REF].

According to PILOTE, simulated drainage along the irrigation period is of 21 mm on SDI120 and 0 for SDI160 in 2009. Simulated yields are still close to average yields as shown in Table 2 which also presents the rainfed treatments One can notice that GY of the rainfed treatments are much lower than the lowest GY of Fig. 7. Although the major part of the root system situated between soil surface and 65cm depth does not profit from the irrigation water, as demonstrated by the sensor analysis, the most remote crops from the pipe still have GY values much higher than which of the rainfed treatments. That reinforces the fact that solely a very limited part of the root system (between 0.7 to 1.2 m depth) contributed to the plant supply. The last column of Table 2 indicates the water productivity calculated as:

WP (kg/m3) = (GY W -GY D )/WAD (5)
In Eq(5), GY W and GY D are grain yields obtained under irrigation and under conditions respectively. Our experimental results seem to be in agreement with the fact that irrigation water is more profitable for roots under SDI120, the latter having a better WP than SDI160. In spite of the high yield variability, this result is compatible with the fact that the lowest lateral spacing insures the best water supply to plants. The quality of the yield simulations from 2007 to 2009 regarding the DI and SDI is noticeable from.Fig. 10 which regroups TDM and GY (at the standard humidity 15%).

Hydrus-2D results

The main objective of this Hydrus-2D application was for drainage and AET evaluations. For that, the simulated cropping cycle was divided into two periods as performed in Mubarak et al., 2009b). The first is from sowing to beginning of July and the second from July to 09/05 in 2008 and to the end of August for 2009, a cropping season much hotter than 2008.

Satisfactory simulations of the soil water content profile  (z) under pipe and under row are presented on Fig. 11 and Fig. 12 for DOY183 and DOY229 of 2009, dates for which  (z) is contrasted (beginning and end of the irrigation season respectively). Simulated drainage is of 20 and 8 mm for first and second period respectively, total drainage representing 11% of the water applied. One can assume these losses could have been reduced if WADs at the beginning of the irrigation season would not have been so high. A comparable application of the model for SDI160 gives much lower drainage losses: 11 mm. Simulated drainage is negligible in 2008 with low WADs for the two SDI treatments. From sowing to DOY = 183 one can assume that no water stress occurred yet according to the soil water status monitoring. The adopted root pattern for that period of 2009 is traced from root density (Mubarak et al., 2009b). At the opposite, that of the end of the cropping cycle, (Fig 12), was adjusted to match less or more with the measured soil water content profile. The resulting rooting pattern clearly showed that a root density notation (in the sense of the Hydrus-2D code) is much higher for the deeper layers and closer to the vertical of the pipe than the one established from our observations just before harvest.

Simulated SWRs on the maximal root depth Px are presented in Fig. 13 for a most conventional case (a pipe at mid distance of two plant rows) as in 2008. For this year, little affected by water stress, the adopted root pattern for Hydrus-2D was based on the observed root density. SWR(Hydrus) is confounded with SWR(PILOTE) just after the beginning of irrigation (end of June). The water amount still available in soil at this period could explain this state of fact, when SWR(PILOTE) should be above measured SWR(under Crop), a hypothesis considerably unquestionably by Fig. 13b having regard to a previous discussion. The perceptible tendency to an over estimation of SWR(under Crop) by Hydrus-2D is more pronounced at the end of the cropping cycle, where the soil water reserve is much more depleted. In spite of the discrepancy between measured and simulated SWR, the decreasing rate is comparable, meaning that simulated and measured AET rate are close. In 2009, SWR(PILOTE) varies between SWR(hydrus) and measured SWR under the drip line for SDI120 (Fig. 14a). SWR(PILOTE) is far from measured SWR under pipe until DOY110 highlighting thus the 2D character of the water transfer process (Fig. 15a) well enough simulated by Hydrus. SWR(Hydrus) is close to measured SWR at the end of the cropping cycle for SDI120 (Fig. 14b) but far for SDI160 (Fig. 15b) where local water stress can be high for such a lateral spacing Fig. 15b shows that SWR(Hydrus) and SWR(PILOTE) for SDI160 are close at the end of the cropping cycle, both overestimating measured SWR, a logical situation for PILOTE but less logical for a 2D soil water transfer such as Hydrus.

Table 2 gives AET values simulated by Hydrus-2D and PILOTE. TDM values are also presented, derived from the empirical equation of Fig. 4. Note that MET = 618 mm for the two treatments of 2008 when MET = 575 mm for SDI120 and 590 mm for SDI160 in 2009. In spite of the root notation fitting at the end of the 2009 season, the AET difference between the two models is still a little higher in 2009 than in 2008 due to higher water stress conditions in 2009 than in 2008 resulting from irrigation strategies and climatic conditions. The fact that MET(2008) is greater than MET(2009) is due to the cropping cycle duration governed by temperature. According to the empirical relationship of Fig. 4, AET underestimations derived from Hydrus-2D result in TDM underestimations for the two contrasted seasons, the discrepancy increasing with the water stress level. These simulation examples seems to point out the impossibility of this Hydrus-2D version to mimic the compensation crop water uptake phenomenon, for simulating AET under severe water stress conditions, water being still less available in the vicinity of roots for SDI160 than for SDI120. Note that a Hydrus-2D simulation for 2009 on the basis of a symmetric root distribution system resulted in a reduction of the AET deviation between the two models of few mm only.

Discussion and Conclusion

The objective of this paper was to show that a 1D crop model such as PILOTE can satisfactory simulate AET and yield under surface and subsurface drip irrigation system whatever the adopted irrigation strategy to reach the target corn yield. Beside this objective, there was the analysis of the role played by the root system through the water uptake process in relation with the soil water availability. For that purpose, the pertinence of a distributed root water uptake model such as which proposed by Hydrus-2D was analysed.

The results demonstrate that even if Hydrus sometimes simulated well enough SWC or SWR thanks to some root density adjustments, the ability of that model for simulating AET under severe water stress conditions is questionable. Such a statement undoubtedly initiated the article of [START_REF] Simunek | Compensated Root Water and Nutrient Uptake Ecological modelling[END_REF] which should be implemented in a future Hydrus version. In the frame of the present work, this statement is strengthened by the narrow relationship existing between AET and yield established for corn on a loamy soil context under a Mediterranean climate. This Hydrus version which cannot account for the compensating root water uptake process would underestimate AET so, at the opposite of a lumped root water uptake model such as which proposed by PILOTE. The latter could give similar actual evapotranspiration rates as those obtained with a compensated water uptake model implemented in Hydrus. But this needs to be verified.

Nevertheless, we have to press on the fact the results refer to a loamy soil where the root system occupies nearly the whole soil domain. Under such conditions, plant can remove its water uptake capabilities to soil regions where water is more easily available. Indeed, the EnviroSMARTt response analysis revealed that lateral extension of humidity does not reach the corn row even for SDI120 in 2008 with a traditional irrigation strategy when a sensitive increase of lateral humidity would be perceptible for layers much deeper than the drip line for irrigation at low frequency in 2009. The relative limited wetting pattern derived from the EnviroSMART response analysis, attest of adaptation capabilities of the root system to restrictive irrigation conditions.

Drainage evaluated by Hydrus-2D, according to the irrigation strategy of 2009 applied to a loamy soil, could have been lower than 10% of the total WAD meaning a better irrigation management.

When restricting the domain of application to a loamy soil, PILOTE can be used for estimating water productivity resulting from irrigation strategies even under drip irrigation. But, at the opposite of Hydrus-2D, PILOTE cannot be recommended for testing another system design on a soil having very different soil properties from a loam. For instance, Hydrus-2D simulations would show that a SDI160 on a sandy loam soil would require a water amount of 440 mm with a necessary frequent irrigation strategy sensitively lower than MET rate on the climatic scenario of 2009. Drainage losses would be high: 150 mm for a grain yield limited to 9.9 T/ha when assuming a possible utilisation of TDM = F(AET) with a harvest index of 0.5. Obviously, these results have to be carefully considered with regards to the fact that we do not know what would be the root system pattern in such a soil type.

Our field experiments carried out on a loamy soil context have shown that adopting a lateral spacing of 1.6 m is a solution which is technically possible even when only 70% of the required water for obtaining the maximal yield would be available. That confirms literature results from USA about the optimal lateral spacing 1.5 m for corn [START_REF] Lam | Subsurface drip irrigation for corn production: a review of 10 years of research in Kansas[END_REF]. But these authors reveal that the SDI system duration should be of 10 to 15 years at least to be more profitable than the pivot system. It is obvious that economical studies would be necessary for justifying the adoption of the SDI system for corn under the SE of France. Nevertheless, the problem of the installation of the root system due to a possible occurrence of dry springs would require a peculiar attention [START_REF] Mailhol | Systèmes d'irrigation et économies d'eau sous climat méditerranéen: Etude expérimentale et modélisation. Actes du séminaire AGDUMED[END_REF]. Utilised in conditions close to our experimental context, an operative model such as PILOTE could help identify economical solutions from simulated irrigation strategies on a climatic series.

Table 1. Soil properties of the loamy soil of Lavalette with the hydraulic parameters used for the Hydrus-2D modelling drawn from Mubarak et al., (2009).s and r are saturated and residual water content,  and n empirical parameters determining the shape of the retention curve, Ks is saturated conductivity, and l a pore connection parameter.

Fig. 2 Fig. 3 .TDMFig. 4 Fig. 5 .Fig. 6 Fig. 7 Fig. 8 .Fig. 9 Fig. 10 .Fig. 11 Fig. 12 Fig. 13 .Fig. 14 .

 234567891011121314 Fig.2 Plan view of the relative position of the neutron probe access tubes, buried pipes, corn rows on SDI120 (a) and SDI160 (b) and location of the EnviroSMART sensors on SDI120 (c) in 2009.

Table 2 .

 2 Grain yield measured and simulated by PILOTE in 2009 with water application depth (WAD) and water productivity (WP).

	Soil Layer (cm)	Clay (%)	Silt (%)	Sand (%)	s	r	(m-1) 	n	(m h -1 ) Ks	l
	0-55	18	42	40	0.36	0	2.73	1.23 1.13 x 10 -02	0.5
	55-90	22	47	31	0.38 0.05 1.30 1.45 5.00 x 10 -03	0.5
	90-150	25	52	18	0.41 0.09 1.90 1.31 2.58 x 10 -03	0.5
	Treatments Measured Mg/ha) Simulated (Mg/ha)	WAD (mm)	WP (kg/m3)
	RF60		7.5		7.0		0		-
			2.5		2.5
	RF75		8.4		8.0		0		-
			3.3		3.1
	SDI160		19.7		20.1		249
			11.8		12.2				3.41
	SDI120		19		19.7		247
			11.8		11.7				3.76
	DI		25		25.6		348
			16		15.8				3.65
	Table 3 Cumulative seasonal AET (mm) and TDM (Mg/ha) evaluated by Hydrus-2D and
	PILOTE from the empirical equation of Fig.4.
	Year Treatment	AET(Hyd) TDM(Hyd) AET(Pil)	TDM(Pil)		TDM(Measured)
	2008 SDI120	545		23.8		585	26.1		26
	2008 SDI160	526		22.8		561	24.7		24
	2009 SDI120	395		15.5		449	18.5		19
	2009 SDI160	435		17.5		480	20.2		20

Author-produced version of the article published in Agricultural WaterManagement, 2011, 98(6), 1033-1044. The original publication is available at http://www.elsevier.com/locate/agwat


Chi-Chung Wang

Tian-Shiang Yang

Phd Candidate

Dynamical Responses of a Valveless Fluid Loop Excited by the Impact of a Compression Actuator

Keywords: valveless pumping, Liebau effect, actuator impact, nonlinear dynamics

Valveless pumping assists in fluid transport in various organisms and engineering systems. In a previous work, to study the actuator impact effects on valveless pumping, we constructed a piecewise-linear lumped-parameter model for a closed-loop system, which consists of two distensible fluid reservoirs connected by two rigid tubes. The preliminary asymptotic and numerical results of that work indicated that the system dynamics is rather rich and complex, and strongly depends upon the driving frequency of the actuator (which periodically compresses one of the distensible reservoirs) and other system parameters.

Here, a more comprehensive numerical study on the aforementioned model is carried out, so as to clarify how many different types of system responses can be excited, and to locate the parameter boundaries within which each type exists.

Moreover, by examining the driving-frequency dependences of a number of characteristic phases (such as the phase of the compression cycle at which the actuator collides with or separates from the distensible reservoir it acts upon), the interrelations between different types of system responses are identified.

INTRODUCTION

In various organisms and engineering systems, the working fluid transport more or less is assisted by valveless pumping. For example, even when the heart's valves fail, blood circulation in the cardiovascular system still is maintained to some extent (see [START_REF] Liebau | Uber Ein Ventilloses Pumpprinzip[END_REF][START_REF] Moser | Impedance Defined flow: Generalisation of William Harvey's Concept of the Circulation -370 Years Later[END_REF]. Various other biomedical systems involving valveless pumping also are discussed in [START_REF] Moser | Impedance Defined flow: Generalisation of William Harvey's Concept of the Circulation -370 Years Later[END_REF]. Meanwhile, in microfluidics, valveless "impedance pumps" have been fabricated by [START_REF] Rinderknecht | A Valveless Micro Impedance Pump Driven by Electromagnetic Actuation[END_REF], and it is expected that they would be particularly suitable for handling sensitive biofluids. One way to achieve valveless pumping is to compress a pliant tube periodically at an asymmetric site from the tube's interfaces to different tubing or reservoirs, so that pressure waves are excited and propagate on the tube. Due to impedance mismatch at the tube ends, the waves would partly reflect, and hence build up a mean pressure head that drives the fluid flow in the tube (see [START_REF] Hickerson | Experimental Study of the Behavior of a Valveless Impedance Pump[END_REF][START_REF] Hickerson | On the Resonance of a Pliant Tube as a Mechanism for Valveless Pumping[END_REF][START_REF] Wen | Design and Characterization of Valveless Impedance Pumps[END_REF].

Despite the existence of a sizable research literature on valveless pumping (for a brief review, see [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF], and references cited therein), a full understanding of its mechanisms has not yet been obtained. To a certain extent, this is because that valveless systems may differ in configuration and operation conditions (e.g., the compression frequency; see [START_REF] Hickerson | Experimental Study of the Behavior of a Valveless Impedance Pump[END_REF]Wen and[START_REF] Wen | Design and Characterization of Valveless Impedance Pumps[END_REF], and therefore rely upon one or more different mechanisms to pump fluids. It is intuitively clear, however, that for a valveless pumping system to produce any mean flow, a necessary condition (but not necessarily a sufficient one; see [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF] is to drive (i.e., compress) the system at an asymmetric site with respect to its configuration. It is also arguable that appropriate theoretical modeling is as important as extensive experimental studies in furnishing a more complete understanding of valveless pumping mechanisms.

As pointed by [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF], in existing mathematical models of valveless pumping, the interaction between the compression actuator and the pliant part of the system in question has often been modeled in rather oversimplified ways. Specifically, in some models the instantaneous deformation of the pliant part of the system is prescribed in a manner not inconsistent (but not absolutely consistent either, as demonstrated by [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF] with the actuator motion (see, e.g., [START_REF] Ottesen | Valveless Pumping in a Fluid-Filled Closed Elastic Tube-System: One-Dimensional Theory with Experimental Validation[END_REF][START_REF] Huang | A Standing Wave Model for Acoustic Pumping Effect in Microchannels[END_REF], while in some Paper Received August, 2010. Revised November, 2010. Accepted December, 2010. Author for Correspondence:Tian-Shiang Yang. other models the temporal variation of the external pressure at the compression site is specified (e.g., Timmermann and Ottesan, 2009). Meanwhile, in the wave-pulse model of [START_REF] Hickerson | On the Resonance of a Pliant Tube as a Mechanism for Valveless Pumping[END_REF], a pair of pressure waves of a prescribed waveform are sent out from the compression site at the driving frequency. While the actuator and the pliant part of the system may indeed interact in a rather simple way under suitable conditions, it was demonstrated by [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF] that there are cases in which the collision of the compression actuator and the system's pliant part is the only cause that produces a mean flow in the system.

As the present work extends that of [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF], let us briefly outline here their model system and major findings. The model system they considered consisted of two distensible reservoirs connected by two rigid tubes, with one of the reservoirs (hereafter referred to as the 0 V reservoir) compressed periodically by an actuator at a prescribed frequency (see Fig. 1). It was assumed that the actuator motion sets the upper bound of the instantaneous volume of the 0 V reservoir. Moreover, during the whole compression cycle, that upper bound is smaller than the free (i.e., unstressed) volume of the 0 V reservoir, and so the nominal duty cycle of the actuator is unity. To suppress wave propagation effects, the driving frequency was assumed to be sufficiently low, and a lumped-parameter model (with constant coefficients) accounting for mass and momentum balance in the system was constructed. Based upon such a model, it was shown that the 0 V reservoir may separate from the compression actuator during some part of the compression cycle if the driving frequency of the actuator exceeds a certain threshold value. The actual duty cycle (i.e., the fraction of the compression cycle during which the actuator actually compresses the system's plaint part) would then become less than unity. In fact, only when that happens would a nonzero mean flow be produced in the fluid loop by the synergetic interaction between configurational asymmetry and the nonlinear effects associated with actuator impact. Moreover, the preliminary asymptotic and numerical results of Yang and Wang (2010) indicated that the system dynamics is rather rich and complex, and strongly depends upon the driving frequency of the actuator and other system parameters. (More technical details of their model and findings will be discussed in the next section.)

In this paper, we shall present the results of a more comprehensive numerical study on the model of [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF]. The purposes are to clarify how many different types of system responses can be excited, and to locate the parameter boundaries within which each type exists. Such a study is important because a good knowledge about the qualitative and quantitative behaviors of the system Fig. 1. Schematic of the model system studied here.

responses is important for designing an efficient valveless pump under certain prescribed specifications. Moreover, by examining the driving-frequency dependences of a number of characteristic phases (such as the particular phase of the compression cycle at which the actuator collides with or separates from the 0 V reservoir), the interrelations between different types of system responses will be identified. Such interrelations clearly are of interest to researchers of nonlinear system dynamics.

MATHEMATICAL MODEL

As in [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF], here we consider a closed-loop valveless pumping system that consists of two distensible fluid reservoirs (such as bellows) connected by two rigid tubes; see Fig. 1. The system is completely filled with an incompressible liquid (say, water) of density ρ and viscosity µ .

The rigid tube on the left-hand side of the system (hereafter referred to as the L-tube) has radius L a and length L l , while that on the right-hand side (the R-tube) has radius R a and length R l . The 0 V reservoir interacts with an external compression actuator, and has a temporally varying volume 0 ( ) V t . Similarly, the instantaneous volume of the other reservoir (the 1 V reservoir) is denoted by 1 ( ) V t .

Constitutive equations

It is assumed here that the physical presence of the actuator sets the maximum admissible volume of the 0

V reservoir at all times. Here we denote such a maximum admissible volume by ( )

a V t .
Accordingly, when 0 ( ) ( ) a V t V t < , the 0 V reservoir will not be in contact with the actuator, and so its volume would vary freely (but in a way that conserves mass and satisfies momentum balance; to be discussed below). For convenience, we shall then say that the system is in the "free mode." On the other hand, whenever the 0 V reservoir is in contact with the actuator-and the system is said to be in the "contact mode"-we shall simply have 0 ( ) ( )

a V t V t = .
As a particular example, it is specified here that ( ) a V t varies sinusoidally, with an angular frequency ω , between min

V ( 0
> ) and max V ( min V > ):
( ) ( )

max min max min 1 1 ( ) cos 2 2 a V t V V V V t ω = + + - . (1)
Moreover, the instantaneous reservoir volumes 0,1 ( ) V t are assumed to vary linearly with the pressure difference across the reservoir walls. Specifically, with 0 ( ) p t being the fluid pressure inside the 0 V reservoir, and ( ) e p t its external pressure (exerted by the actuator), one has

{ } 0 0 0 0 ( ) ( ) ( ) f e V t V C p t p t = + ⋅ - , (2)
where 0 f V is the free (i.e., unstressed) volume of the 0 V reservoir, and 0 C is a (constant) compliance coefficient. It is assumed here that the actuator can only exert a positive pressure on the 0 V reservoir, so that 0 e p > when the system is in the contact mode, and 0 e p = in the free mode. Meanwhile, since the 1 V reservoir has an internal pressure 1 ( ) p t but no external pressure, its constitutive relation is written as

1 1 1 1 ( ) ( ) f V t V C p t = + ⋅ , (3)
with 1 C being a constant compliance coefficient.

Mass and momentum balances

Now, denoting the volumetric flowrates (from the 0 V reservoir to the 1 V reservoir) in the L-and R-tubes by ( ) , 4b) where the overdots represent derivatives of the state variables with respect to time t . Meanwhile, to suppress all nonlinearities except that arising from the interactions between the actuator and the 0 V reservoir, the momentum equations for the L-and R-tubes are given by T I R = . Similarly, by Eq. (5b), the intrinsic time scale for the variation of ( )

L Q t and ( ) R Q t , respectively, mass conservation requires that 0 L R V Q Q = -- & , 1 L R V Q Q = + & , ( 4a
0 1 L L L L I Q R Q p p + = - & , ( 5a
) 0 1 R R R R I Q R Q p p + = - & , (5b) where
R Q t when the R-tube stands alone is / R R R T I R = .

Initial conditions

As in [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF], here we shall consider the particular case that both the 0 V and 1 V reservoirs are filled to their free volumes before the compression actuator is installed. The total fluid volume in the two reservoirs,

0 1 0 1 ( ) ( ) tot f f V V t V t V V = + = + , therefore is invariant
with time as can be readily deduced from Eqs. (4a) and (4b). Suppose then that an actuator is installed and compresses the 0 V reservoir to a certain extent before it starts moving.

In terms of the "pre-compression ratio" ( )

0 max 0 / f f V V V γ = - ( 0 > ),
the initial volumes of the 0 V and 1 V reservoirs therefore are given by 0 0

(0) (1 ) f V V γ = - and 1 1 0 (0) f f V V V γ = +
, respectively. Also, Eq. ( 1) can be rewritten as

( ) 0 ( ) (1 ) 1 cos a f a V t V V t γ ω = --∆ - , (6)
where max min

( ) / 2 a V V V ∆ = - . Meanwhile, one has (0) 0 (0) L R Q Q = =
, and it can be calculated from Eqs.

(2) and (3) that initially the uniform internal pressure of the system is 0

0 1 1 (0) / (0) f p V C p γ = =
, and the external pressure exerted by the actuator on the 0 V reservoir is

0 01 (0) / e f p V C γ = , with the "characteristic compliance" of the system 01 C given by 01 0 1 1/ 1/ 1/ C C C = + .
With the initial conditions of the system specified above, Eqs. ( 2)-( 6) can be integrated numerically (say, by use of the fourth-order Runge-Kutta method). However, as shown in [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF], for such a pre-compressed system if the driving frequency ω is below a well defined threshold value th ω , the 0 V reservoir would never separate from the compression actuator, and the perfectly linear governing equations would produce zero mean flowrate in the fluid loop. As a matter of fact, since the purpose of this work is to examine the effects of actuator impact upon valveless pumping, the pre-compressed system is deliberately chosen, such that a nonzero mean flow in the system can only be produced by the interactions between the actuator and the 0 V reservoir.

Figures 2(a) and 2(b) summarize the possible interactions between the actuator and the 0 V reservoir.

Since such interactions have been discussed in detail in [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF], we shall only briefly highlight their physical meanings below.

Mode decision rules

First, suppose that at the present time step the system is in the contact mode. As explained above, one will then have 0 ( ) ( )

a V t V t =
, and it follows from

Eq. ( 4a) that ( )

( ) ( ) R a L Q t V t Q t = - - & . Moreover, the pressure difference 0 1 ( ) ( ) p t p t - can be eliminated from Eqs. (5a) and (5b), yielding L L L L I Q R Q + = & ( ) R R R R R a L I Q R Q I V Q + = - + & & && ( ) R a L R V Q - + & . With
the particular choice of ( ) a V t given by Eq. ( 6), this result can be written as 2 ( ) ( )

( sin cos ), L R L L R L a R R I I Q R R Q V R t I t ω ω ω ω + + + = ∆ + & (7
) from which the value of L Q at the next time step can be readily calculated by use of the fourth-order Runge-Kutta method. (Incidentally, Eq. ( 7) also shows that in the contact mode the characteristic time scale for the flowrate variation in the closed fluid loop is ( ) /( )

LR L R L R T I I R R = + +
.) The values of other state variables at the next time step then can be calculated in turn following the procedures summarized in Fig. 2(a). If a negative value of the external pressure e p is obtained at the next time step, it is interpreted that the system switches to the free mode (otherwise the system stays in the contact mode at the next time step). Of course, when a negative value of e p arises, one would also need to determine the precise "separation time" ( s T , when 0 e p = ) by interpolation; see Fig. 2(a). Now, if the system is in the free mode at the present time step, then the external pressure ( ) 0 e p t = . Accordingly, it is deduced from Eqs. ( 2) and ( 3) that 0

1 0 0 01 ( ) ( ) [ ( )] / f p t p t V V t C - = - - , which
can be substituted into Eqs. ( 5a) and (5b) to yield It should also be noted that if the system currently is in the free mode, there are more decisions to be made regarding the system state at the next time step. Specifically, in the free mode, the volume of the 0

0 0 01 [ ( )] / L L L L f I Q R Q V V t C + = - - & , ( 8a
) 0 0 01 [ ( )] / R R R R f I Q R Q V V t C + = - - & . ( 8b
V reservoir may reduce to zero or increase to meet its instantaneous upper bound ( ) a V t set by the actuator motion; the latter is interpreted as a collision of the 0 V reservoir with the actuator. In the computations, if 0 ( ) V t becomes negative in the next time step, we shall determine first the true "depletion time" 0 T , with 0 0 ( ) 0 V T = , by interpolation. Moreover, it is assumed that as soon as the 0 V reservoir is completely depleted, the fluid loop is blocked, so that 0 0

( ) 0 ( ) L R Q T Q T = =
and, of course,

1 0 0 1 ( ) f f V T V V = + . For 0 t T + =
, the system is restarted with the state variable values given above, and then continues its journey in the free mode (see Fig. 2(b)).

Consider next the possibility of collision between the 0 V reservoir and the actuator. Again, in the computations it may occur that at the next time step 0 ( ) ( )

a V t V t > ,
L L c R R c a c L c L R I Q T I Q T V T Q T I I - - + - + = + & , ( 9a
) ( ) [ ( ) ( )] ( ) R R c L L c a c R c L R I Q T I Q T V T Q T I I - - + - + = + & . ( 9b
)
Another important calculation here is the system mode following the collision. Specifically, after the collision, the system may switch to the contact mode, or stay in the free mode. To decide what would happen, in the computations we assume first that the system stays in the free mode at

( ) ( ) a V t V t <
, we shall take that the system really would stay free; otherwise the system would switch to the contact mode (see Fig. 2(b)).

Clearly, the overall mathematical model involves rather intricate decision rules (summarized in Figs. 2(a) and 2(b)), and its solution generally has to be obtained numerically. However, for a set of realistic sample parameter values (the same as that specified in the next section), [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF] identified a small dimensionless compliance parameter, and therefore were able to carry out a series of asymptotic calculations to obtain accurate estimates for some critical driving frequencies across which the system responses undergo qualitative changes. Furthermore, numerical solutions of the model for three representative cases (with min max / V V = 0.6, 0.33, and 0.15) were carefully discussed to elucidate the complex dynamics of the model system, and to illustrate the frequency dependence of the mean flowrate and mean pressures in the system. Their numerical results indicated that there may exist steady periodic system responses having a period that is a multiple of the compression cycle period at higher driving frequencies. In fact, under suitable conditions, seemingly chaotic system responses also arose in their computations. However, the transition routes of the various types of system responses were not fully discussed by [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF]. Moreover, for steady periodic system responses repeating themselves in one compression cycle, [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF] had not determined the precise boundaries between different types of system responses. These, then, are the main purposes of the present work.

SAMPLE PARAMETER VALUES

In our computations, a sample system is envisioned which has realistic dimensions and material properties. Specifically, it is supposed that the 0

V and 1 V reservoirs are made of bellows.

The free (i.e., unstressed) radius and height of the 0 V reservoir are 0 25 r = mm and 0 33.3 h = mm, respectively, while the corresponding dimensions of the 1

V reservoir are 1 50 r = mm and 1 66.7 h = mm. The free volumes of the two reservoirs therefore are

L R LR L R C R R T R R ε = + 5 1.03 10 1 - = ×
, and therefore were able to carry out some asymptotic analyses. Note also that, as in [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF], here the maximum value of ( )

a V t is taken to be max 0 0.6 f V V = .
(The numerical results with other choices of max V would appear to be qualitatively similar to that discussed in the next section, and therefore will not be presented here.) Meanwhile, to avoid degenerate limiting cases, the minimum value of ( )

a V t , min V ,

RESULTS AND DISCUSSION

Classification of system response

Let us now discuss the numerical results of our impact dynamics model for a model valveless pumping system having the parameter values specified above. First, Fig. 3 shows the parameter boundaries for various types of system responses observed in our computations. Note that, to bring our some dynamical similarities, the numerically calculated critical frequencies are normalized by Basically, in the parameter range shown in Fig. 3, when the driving frequency ω of the actuator does not exceed the "stability-margin frequency" st ω , five different types of system responses are observed, namely the "pure contact" (labeled by abbreviation PC in Fig. 3), "contact-free" (CF), "purely free" (PF), "contact-free with 0 V reservoir depletion" (CF0),

and "purely free with 0 V reservoir depletion" (PF0) responses.

The PC responses exist when ω is below the threshold frequency ( th ω ), in which the 0 V reservoir never separates from the actuator and no mean flowrate is produced in the fluid loop. When the volume ratio min max / V V is larger, the CF responses arise when ω exceeds th ω but is lower than the onset frequency PF ω of the PF responses. In the CF responses, the system is in the contact mode for a certain portion of the compression cycle, and is in the free mode for the rest of the compression cycle; see For smaller values of the min max / V V ratio, when ω exceeds th ω , the 0 V reservoir also would separate from the actuator during some portion of the compression cycle, so that the CF responses would arise. However, beyond a certain critical driving frequency 0 CF ω , the 0 V reservoir would be completely depleted at a certain phase in the free-mode portion of the CF responses. Such responses (with complete 0 V reservoir depletion) are referred to as the CF0 responses; see Fig. 4(c). seemingly chaotic at the stability-margin margin frequency st ω . Also, since the system essentially is restarted at the instant when the 0 V reservoir is completely depleted (i.e., when 0 ( ) 0 V t = ), and therefore intuitively would be somewhat more stable, the stability-margin frequency for min max / V V < 0.3342 is much higher than that for min max / 0.3342 V V > ; see Fig. 3(a).

Frequency dependences of various characteristic phases

To begin with, let us consider the case of θ at which the actuator collides with the 0 V reservoir.

To examine the driving-frequency dependences of such characteristic phases, and bring out their interrelations at the same time, we find it convenient to use polar plots.

Specifically, in Fig. 5

th f f - = Hz)
. It is also interesting to note that, since the system responses undergo a qualitative change (from PF to PF0) at such a critical frequency, the r θ curve then makes a sharply turn there; see Fig. 5(b). Meanwhile, since the system essentially is restarted after the 0 V reservoir is completely depleted, the stability-margin frequency st f is significantly increased (compared with that in the previous case of min max / 0.6 V V = ). Still, the PF0 responses would become seemingly chaotic (which will be discussed later) when f exceeds again, the PF0 responses would eventually become seemingly chaotic when f exceeds st f .

Finally, for the case of min max / 0.15 V V = , the characteristic phases are plotted against the driving frequency f in Fig. 5(d). (Note that, as the range of the driving frequency f is now much larger than that in previous cases, here the value of f is plotted in logarithmic scales.) Compared with the previous case (with min max / 0.325

V V =
), here the transition of the system responses are much simpler, namely from PC to CF (at

th f f = ), CF0 (at 0 CF f f = ), PF0 (at 0 PF f f =
) in turn, and finally to the seemingly chaotic responses (at st f f = ). Also, since all the important features of the characteristic phase curves have been explained above, they are readily understandable in Fig. 5(d), and therefore will not be over-discussed here.

Seemingly chaotic system responses

We now briefly discuss some interesting features of the seemingly chaotic responses of the system when the driving frequency f exceeds the stability-margin frequency st f . Note that, as such responses seem to be chaotic, they are quite sensitive to the specific parameter values of the system. However, for the purpose of demonstrating a number of qualitative features of such responses, here we shall only consider the case of min max / 0.33 V V = . It is also important to note that, for such responses the temporal means (over one compression cycle) of the flowrates in the system would vary from cycle to cycle. So, to extract some orders out of the seemingly chaotic responses, for each parameter setting the numerical computations are run for 2000 compression cycles, and the temporal means of the aforementioned system state variables are computed for each compression cycle. Of course, during the first few cycles, the variations in such temporal means may result partly from the initial transients. However, we found that the initial transient effects practically are negligible after a few dozen cycles. Also, it suffices to keep track of the last 50 cycles (of the 2000 cycles of computation) to bring out the general pictures of the system responses.

To illustrate the frequency dependences of the system responses for As can be readily seen, for the "2T" and "3T" responses, the mean flowrates of the 50 cycles plotted in Fig. 6(b) appear to repeat every two cycles and every three cycles, respectively, and that is exactly why the responses are named as such. In fact, while being somewhat more difficult to observe in Fig. 6(a), some more regularities in the system responses can be observed in Fig. 6 Of course, there certainly are other less obvious regularities of the system responses in the frequency range plotted in Figs. 6(a) and 6(b), and much much more at higher driving frequencies. Meanwhile, outside the particular frequency windows within which such regular responses exist, the seemingly chaotic responses can get as complicated as one can imagine. However, as a detailed analysis of such response clearly is non-trivial, and is beyond the scope of this work, it will not be further pursued. Instead, let us point out an additional observation in Fig. 6(a). Specifically, it is seen in Fig. 6(a) that the V 0 reservoir may stay in contact with the compression actuator only when the driving frequency is lower than about 135 Hz. Above that frequency, the separation-collision interaction between the V 0 reservoir and the actuator would no longer be possible (and the separation and collision phases cease to appear in Fig. 6 (a)). What happen then are that, whenever the V 0 reservoir collides with the actuator they separate from each other immediately afterwards, and that each collision of them would be so strong to completely deplete the V 0 reservoir. As a result, in Fig. 6(a), for f > 135 Hz, only the rebound and zero volume phases exist and are plotted there.

Finally, to conclude this discussion, let us speculate a few possible applications of the system responses with driving frequencies [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF]; the interested reader is referred to that work for further details. Note, however, that since the mean flowrates for st f f > vary from cycle to cycle, and that their overall (i.e., long-term) averages actually have the same order of magnitude as that for st f f < . This finding may have an important engineering bearing. Specifically, if one would like to design a micro-pump with a stable flowrate output, it is advisable to choose a driving frequency f that is lower than the stability-margin frequency st f , because increasing f would likely induce the seemingly chaotic system responses but would not increase the long-term average of the system flowrate. However, if one wishes to design a micro-mixer, the large yet irregular flowrate oscillation for st f f > may be useful for enhancing the mixing effects. It is important to note, however, that some further works still are needed to assess the feasibility of such a possibility.

CONCLUDING REMARKS

In this work, using the lumped-parameter model of [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF] for a valveless pumping system with actuator impact effects, we have carried out a more comprehensive numerical study than that discussed in [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF]. For periodic system responses repeating themselves in one compression cycle, five different types of responses are found, which are referred to above as the pure contact (PC), contact-free (CF), purely free (PF), contact-free with 0

V reservoir depletion (CF0), and purely free with 0 V reservoir depletion (PF0) responses (and illustrated in Fig. 4). By examining the dependences on the volume ratio min max / V V of a few critical frequencies at which the system responses undergo a transition in type, the parameter boundaries within which such responses exist have been clarified (see Fig. 3). Moreover, for a few particular values of min max / V V , the driving-frequency dependences of a number of characteristic phases, namely the separation phase s θ , collision phase c θ , rebound phase r θ , and 0

V reservoir depletion phase 0 θ (see Fig. 5), the interrelations of the different types of system responses have been identified. It is believed that these results not only are useful for understanding the nonlinear impact dynamics of valveless pumping systems, but also would help optimize the design of practical valveless pumps with certain prescribed specifications. It has also been demonstrated in this work that, for driving frequencies f higher than the stability -margin frequency st f , the system responses generally do not repeat themselves after one compression cycle. Instead, they appear to be "chaotic."

For such seemingly chaotic system responses, several interesting features have been observed. For example, we have found several driving-frequency windows within which the system responses appear to repeat themselves after an integral multiple of the compression period (see Figs. 6 and7). Moreover, the mean flowrates (in each cycle) for st f f > generally are greater than that for st f f < by several orders of magnitude, while their overall (i.e., long-term) average has the same order of magnitude as that for st f f < . This finding suggests that, if one would like to design a valveless pump (say, for certain microfluidic or biomedical applications) with a stable flowrate output, it would be advisable to choose a driving frequency f that is lower than the stability-margin frequency st f . However, if one wishes to design a micro-mixer instead, the large yet irregular flowrate oscillation for st f f > may be useful for enhancing the mixing effects. And of course, some further works still are needed to assess the feasibility of such possibilities.

Finally, to conclude this paper, we would also like to point out that the impact dynamics model used here (and in [START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF] has suppressed many other mechanisms that may produce significant valveless pumping effects (such as that of impedance pumps; see, for example, [START_REF] Hickerson | On the Resonance of a Pliant Tube as a Mechanism for Valveless Pumping[END_REF][START_REF] Hickerson | Experimental Study of the Behavior of a Valveless Impedance Pump[END_REF][START_REF] Huang | A Standing Wave Model for Acoustic Pumping Effect in Microchannels[END_REF][START_REF] Wen | Design and Characterization of Valveless Impedance Pumps[END_REF]. The model has been designed in this particular way to demonstrate the complicated actuator impact dynamics of valveless pumping systems with the highest possible clarity. Given the many interesting findings discussed above, in future work it certainly is worth pursuing to design a suitable experimental apparatus for their verification. Meanwhile, in many realistic microfluidic or biomedical applications, one may also need to take into account the synergetic interactions of the actuator impact effect with other valveless pumping mechanisms, which is also worth pursuing in future work.

NOMENCLATURE

  the next time step can be calculated by integrating Eqs. (4a), (8a), and (8b). The values of other state variables at the next time step then can be calculated in turn following the procedures summarized in Fig. 2(b).

  Fig. 2. Numerical solution procedures and mode decision rules.

  and interpolation is needed to determine the true collision time c flowrates in the rigid tubes to vary discontinuously. Specifically, as calculated by[START_REF] Yang | Effects of Actuator Impact on the Nonlinear Dynamics of a Valveless Pumping System[END_REF], one would have

  If it is calculated at the tentative time step that 0

  the two fluid reservoirs (bellows) are reinforced by elastic springs so that they have the same effective spring constant of , is taken to be the working fluid in the system. Using then the Poiseulle's law, the flow resistance coefficients, L R R appearing in Eqs. (5a)and (5b) are related to the corresponding tube lengths , 5a) and (5b) are taken to be proportional to the mass of the fluid in the L -and R

  while their dimensional values also are shown in Fig.3to give some ideas of their actual magnitudes).

  Fig. 4(a) for an example. Meanwhile, in the PF types of system responses and the critical frequencies determining the parameter boundaries of such responses are detailed in the text.responses, the 0 V reservoir collides with the actuator at a particular phase of the compression cycle and separates from the actuator right afterwards; see Fig.4(b). If ω is further increased beyond the stability-margin frequency st ω , the system responses then become much more complicated, and will be discussed later.

  turn, would become the PF responses at PF ω ω = . Also, when ω exceeds st ω , more complicated system responses would be excited. So, for

  variable is simply the phase of the actuator motion. It is clearly seen in Fig.5(a) that the s θ and c θ curves both start at the threshold frequency th f f = , with the common value of about 180 o . This, in fact, is easy to understand, because for th f f < the 0 V reservoir would never separate from the actuator (and therefore be impossible to collide with it). As the driving frequency (i.e., the radial variable) increases, the value of s θ decreases, while that of c θ increases. Accordingly, the actual duty cycle (i.e., the fraction of the compression cycle in which the actuator is in solid contact with the 0 V reservoir) decreases with increasing f . in Fig.5(a)), the actual duty cycle becomes zero, in which case the 0 V reservoir collides with the actuator at a certain phase of the compression cycle and rebounds from the actuator right afterwards. So, logically, in such cases the separation s θ also is the collision phase c θ . As a result, for f ∈ responses would not repeat themselves in one compression cycle, and would appear to be chaotic. Such responses will be discussed later. Consider now the case of min max / that in Fig.5(a). However, in the present case, when f exceeds PF f , the system responses would first experience a complete depletion

  Frequency dependences of various characteristic phases of the system responses for driving frequencies

  remains useful to trace the variations of the critical phase angels s θ , c θ , r θ , and 0 θ , as shown in Fig. 6(a). Note, however, that such phases generally vary from cycle to cycle, and that the results for the last 50 cycles are plotted. So, for a given driving frequency f (which can be represented by a circle in Fig. 6(a)), a certain critical phase (say r θ ) may have multiple values and hence be represented by several data points. Moreover, for st f f > , the critical phases are extremely sensitive to a small variation in f . As a result, the data points in Fig. 6(a) appear to be scattering, and generally do not form continuous curves (while they do in Figs. 5(a)-5(d)). (Note also that here the computations are carried out at the driving-frequency increment of 1 Hz.) However, some characteristic phases do seem to form continuous curves within a number of driving-frequency windows. For example, in the frequency window bounded approximately by 109 Hz and 177 Hz (labeled with "2T" in Fig. 6(a)), the separation phase ( s θ ) and collision phase ( c θ ) appear to form two continuous curves, which start at about 109 f = Hz and merge into a single curve of rebound phase ( r θ ) at about f = 135 Hz. Meanwhile, the zero volume phase ( 0 θ ) also appears to form a continuous curve throughout this frequency window. Moreover, in the frequency window bounded approximately by 226 Hz and 266 Hz (labeled with "3T"), the rebound phase ( r θ ) and zero volume phase

( 0 θ

 0 ) appear to form two continuous curves. Some insights about the dynamical behaviors of the "2T" and "3T" responses can be gained by plotting the mean flowrate in the L-tube, various characteristic phases. (a) and mean flow rate in the L-tube (b) for system responses with driving frequencies (subtracted by the threshold frequency th f ); see Fig. 6(b). Two examples of such responses are indicated by lines (a) and (b) in Fig. 6(b).

  (b); two examples are indicated by lines (c) and (d) there. Specifically, on lines (c) and (d), the mean flowrates of the 50 of (a) 150 Hz, (b) 250 Hz, (c) 90 Hz, and (d) 205 Hz, respectively. plotted cycles appear to repeat every three cycles and every five cycles, respectively. The regularities of such responses also can be seen more clearly by plotting them in time domain; see Figs. 7(a) -(d).

  Fig. 6(b). But the results for

ACKNOWLEDGEMENT

The authors gratefully acknowledge the Taiwan National Science Council for supporting this work through grant NSC96-2628-E-006-241-MY3. They would also like to thank Professors M.-S. Ju, S.-Y. Lee, K.-S. Chen, and C.-Y. Wen of NCKU for a number of fruitful discussions on this work and other related topics.

).

Moreover, when ω exceeds a even higher critical frequency 0 PF ω , the system response will become purely-free with complete depletion of the 0 V reservoir at a certain phase of the compression cycle. Such responses are refereed to as the PF0 responses; see Fig. 4(d). And again, when ω is further increased beyond the stability-margin frequency st ω , the system responses become more complicated (to be discussed later).

It is also important to note that, for intermediate values of the min max / V V ratio, the transition of system responses are somewhat more complicated. Specifically, as shown in Fig. 3 As the driving frequency of the system is further increased, the CF0 responses become the PF0 responses, which then lose their stability at st ω ω = . In order to clarify the interrelations of the various types of system responses, in the following subsection we shall pick four particular values of min max / V V , namely 0.6, 0.33, 0.325, and 0.15 (marked by the Roman numerals I, II, III, and IV, respectively, in Fig. 3), and carefully examine the frequency dependences of a few characteristic phases of the responses (such as the particular phase at which the 0 V reservoir separates from or collides with the compression actuator).

Before we do that, however, an additional feature of Fig. 3