domain: info.info-au math spi sdv

text:
Ludovic Mailleret
email: ludovic.mailleret@sophia.inra.fr

Frédéric Grognard
email: frederic.grognard@inria.fr

Global stability and optimisation of a general impulsive biological control model

Keywords: Predator-prey dynamics, impulsive model, global stability, optimisation

An impulsive model of augmentative biological control consisting of a general continuous predator-prey model in ordinary differential equations, i.e. a meta-model, augmented by a discrete part describing periodic introductions of predators is considered. The existence of an invariant periodic solution that corresponds to prey eradication is shown and a condition ensuring its global asymptotic stability is given. An optimisation problem related to the preemptive use of augmentative biological control is then considered. It is assumed that the per time unit budget of biological control (i.e. the number of predators to be released) is fixed and the best deployment of this budget is sought in terms of release frequency. The cost function to be minimised is the time needed to reduce an unforeseen prey (pest) invasion occuring at a worst time instant under some harmless level. The analysis shows that the optimisation problem admits a countable infinite number of solutions. An argumentation considering the required robustness of the optimisation result with respect to the invasive prey population level and to the model parameters is then conducted. It is shown that the cost function is decreasing in the predator release frequency so that the best deployment of the biocontrol agents is to carry out as frequent introductions as possible.

Introduction

Biological control aims to reduce the populations of harmful organisms by utilising their natural enemies, the biological control agents. It can be used to control vectors in vector-borne diseases like malaria or invasive species (may they be animals or plants), but it is mostly dedicated to the control of insect pests, especially at the agricultural cropping system scale [START_REF] Murdoch | Biological control in theory and practice[END_REF]. Four main biological control strategies can be identified (see e.g. [START_REF] Eilenberg | Suggestions for unifying the terminology in biological control[END_REF][START_REF] Roderick | Genes in new environments: Genetics and evolution in biological control[END_REF] and references therein):

• conservation biological control: biocontrol agents are already present in the system and are favoured through various means (habitat management, cropping practices...).

• classical biological control: exotic biocontrol agents are introduced once in the system and expected to establish and persist indefinitely to provide long term pest control.

• inoculation biological control: biocontrol agents are introduced (once or repeatedly) and are expected to reproduce to provide pest control for an extended period, but not permanently.

• inundative biological control: biocontrol agents are repeatedly introduced in order to eradicate the pest population. Pest control is expected to be achieved by the released organisms themselves and not by their potential progeny.

As noted by Eilenberg et al. [START_REF] Eilenberg | Suggestions for unifying the terminology in biological control[END_REF], it might be difficult to distinguish between repeated inoculation biological control and inundative biological control. Indeed, biocontrol agents used in an inundative program are very likely to reproduce and the distinction between the "inundative" effects and the "inoculative" effects may not be that clear.

Hence, "augmentative" biological control is commonly used to denominate as well repeated inoculation biocontrol or inundative biocontrol, though it is not the perfect denomination either [START_REF] Eilenberg | Suggestions for unifying the terminology in biological control[END_REF].

In this contribution we focus on this "augmentative" method, which is fairly well modelled by an impulsive system of ordinary differential equations (a continuous system affected by sudden state modifications at some discrete time moments). The intrinsic dynamics of the two interacting populations (pests and biocontrol agents) are described by the continuous part of the system and the repeated (periodic) biocontrol agents introductions by the discrete one.

Early examples of the use of impulsive models in biology are [START_REF] Ebenhoh | Coexistence of an unlimited number of algal species in a model system[END_REF][START_REF] Funasaki | Invasion and chaos in a periodically pulsed mass-action chemostat[END_REF] in chemostat modelling. Since then, they have received considerable attention in the area of epidemiology [START_REF] Onofrio | Stability properties of pulse vaccination strategy in SEIR epidemic model[END_REF][START_REF] Hui | Impulsive vaccination of SIR epidemic models with nonlinear incidence rates[END_REF][START_REF] Shulgin | Pulse vaccination strategy in the SIR epidemic model[END_REF], cancer treatment [START_REF] Onofrio | Tumour eradication by antiangiogenic therapy: analysis and extensions of the model by hahnfeldt et al[END_REF][START_REF] Onofrio | A general framework for modeling tumor-immune system competition and immunotherapy: Mathematical analysis and biomedical inferences[END_REF][START_REF] Lakmeche | Bifurcation of non trivial periodic solutions of impulsive differential equations arising chemotherapeutic treatment[END_REF], populations dynamics [START_REF] Liu | Recent progress on stage-structured population dynamics[END_REF][START_REF] Liu | Complex dynamics of holling type ii lotka-volterra predator-prey system with impulsive perturbations on the predator[END_REF][START_REF] Tang | Multiple attractors in stage-structured population models with birth pulses[END_REF][START_REF] Liu | The dynamics of a prey-dependent consumption model concerning impulsive control strategy[END_REF][START_REF] Negi | Dynamics in a Beddington-DeAngelis prey-predator system with impulsive harvesting[END_REF] or integrated pest management [START_REF] Liu | The dynamical behaviors of a Lotka-Volterra predator-prey model concerning integrated pest management[END_REF][START_REF] Nundloll | The effects of partial crop harvest on biological pest control[END_REF][START_REF] Nundloll | Impulsive biological control with intra-predatory interference: global stability and convergence[END_REF][START_REF] Liu | Dynamic complexities of a Holling I predator-prey model concerning periodic biological and chemical control[END_REF][START_REF] Tang | Models for integrated pest control and their biological implications[END_REF]; we refer to e.g. [START_REF] Mailleret | A note on semi-discrete modelling in the life sciences[END_REF] for a recent overview. Though there are numerous studies on bio-inspired impulsive models, to our knowledge only a few of these are concerned with optimisation: see [START_REF] Dong | Optimal harvesting policies for periodic gompertz systems[END_REF][START_REF] Erdlenbruch | Cyclical versus non-cyclical harvesting policies in renewable resource economics[END_REF][START_REF] Meng | A new mathematical model for optimal control strategies of integrated pest management[END_REF][START_REF] Tang | Optimal impulsive harvesting on non-autonomous beverton-holt difference equations[END_REF][START_REF] Xiao | Optimal impulsive control in periodic ecosystem[END_REF]. Here we focus both on a global convergence result for a rather general predator-prey model with periodic impulsive predator introductions and on a problem of optimisation on how to best deploy these introductions in the context of crop protection through biological control. Our results follow and generalise earlier ones that were obtained in a more restricted context [START_REF] Mailleret | Optimal release policy for prophylactic biological control[END_REF].

In this paper, we first present the impulsive predator-prey model of augmentative biological control that is considered. The hypotheses assumed on the biological part (continuous system) are qualitative so that they can encompass a large number of classical biological/ecological functions. The existence of an invariant periodic prey-eradication solution is proved and a sufficient condition that ensures its global asymptotic stability is given.

Then, we focus on an optimisation problem in the perspective of preemptive/preventative use of augmentative biological control. This departs from most classical biological control procedures which are of a more "feedback" type: biocontrol agents are released once the pests have been detected only. Such a preventive approach is also important because it appears to achieve more acceptable pest control especially for high valued agricultural crops like e.g. orchards, vegetables, mushrooms or ornamentals, that are very sensitive to the slightest pest outbreak (see [START_REF] Fenton | Evaluating the efficacy of entomopathogenic nematodes for the biological control of crop pests: A nonequilibrium approach[END_REF][START_REF] Skirvin | Modelling the effects of plant species on biocontrol effectiveness in ornamental nursery crops[END_REF][START_REF] Fenton | Optimal application strategies for entomopathogenic nematodes: integrating theoretical and empirical approaches[END_REF] for theoretical/simulation studies and [START_REF] Ehler | Potential for augmentative biological control of black bean aphid in california sugarbeet[END_REF][START_REF] Jacobson | Suppressing establishment of frankliniella occidentalis pergande (thysanoptera : Thripidae) in cucumber crops by prophylactic release of amblyseius cucumeris oudemans (acarina : Phytoseiidae)[END_REF][START_REF] Jacobson | Compatibility of beauveria bassiana (balsamo) vuillemin with amblyseius cucumeris oudemans (acarina : Phytoseiidae) to control frankliniella occidentalis pergande (thysanoptera : Thripidae) on cucumber plants[END_REF][START_REF] De | Biological control of thrips on ornamental crops: Interactions between the predatory mite neoseiulus cucumeris (acari : Phytoseiidae) and western flower thrips, frankliniella occidentalis (thysanoptera : Thripidae), on cyclamen[END_REF][START_REF] Hulshof | Life history characteristics of frankliniella occidentalis on cucumber leaves with and without supplemental food[END_REF] for real life experiments and a biological perspective).

It is assumed that a fixed budget ensuring pest eradication is allowed and that only its deployment (a tradeoff between frequent-small or rare-large releases) may be modified. The cost function to be minimised is the time taken to reduce 2 an unforeseen pest invasion under some harmless level and the input is the impulsive release frequency.

It is shown that the local optimisation problem has a countable infinite number of solutions and, in a second step using an argumentation based on the needed robustness of our result, it is proved that the locally optimal solution is to choose as frequent (and thus as small) releases as possible. It is also shown that this strategy gives a sub-optimal solution to the global optimisation problem. Moreover, the higher the frequency, the lower the cost function, so that even if one cannot achieve a very high frequency, it is always interesting to choose the highest possible one. A Monte Carlo like numerical simulation illustrates our analytical optimisation result and the article is concluded by some recommendations to biological control practitioners that follow from our analysis.

Model Description and Analysis

A general impulsive biological control model

In this paper, we assume that the tri-trophic system (crop -pests -biocontrol agents) may be well represented by its bi-trophic simplification representing the prey (pest) predator (biocontrol agent) interactions. The underlying assumption is that the prey population remains at moderate levels so that the crop is not limiting for prey growth and thus do not need to be taken into account in the model. We consider the simple case in which no intra specific interactions within the predator population affects the predator prey interaction, though some in the prey population may occur (see Remark 1 below). According to classical predator-prey modelling, we get the following two-dimensional

model    ẋ = f (x) -g(x)y ẏ = h(x)y -my (1)
with x denoting the prey density and y the predator density. f (x) denotes the prey's growth rate, g(x), h(x) and m denotes the predators' functional response, numerical response and mortality rate, respectively. Since biological processes are always difficult to model, we only assume general qualitative hypotheses on the functions f (.), g(.) and h(.) (see Hypothesis 1 below). Thus, we consider a family of models, i.e. a "meta-model" sensu [START_REF] Onofrio | Tumorimmune system interaction: Modeling the tumor-stimulated proliferation of effectors and immunotherapy[END_REF][START_REF] Onofrio | Metamodeling tumorimmune system interaction, tumor evasion and immunotherapy[END_REF], rather than a specific mathematical model. Such an approach guarantees that our analyses apply to a wide range of biological systems, may they be directly related to biological control or not. Indeed, many are the biological systems that fall in the metamodelling [START_REF] De | Biological control of thrips on ornamental crops: Interactions between the predatory mite neoseiulus cucumeris (acari : Phytoseiidae) and western flower thrips, frankliniella occidentalis (thysanoptera : Thripidae), on cyclamen[END_REF]: examples may for instance be found in immunology, epidemiology or cancer modelling.

Hypotheses 1 (H1): Let f (.), g(.) and h(.) be locally Lipschitz functions on R + such that:

(i) f (0) = 0 (ii) g(0) = 0, g (0) > 0 and ∀x > 0, g(x) > 0 (iii) the function f (x) g(x)
is upper bounded (for x > 0)

(iv) h(0) = 0 and ∀x > 0, h(x) ≥ 0 Remark 1: A large part of the predator-prey functions encountered in the literature fit these hypotheses: (H1-i) only indicates that no spontaneous generation of pests (prey) is possible; f (.) may be constant, linear, logistic or model an Allee effect etc...(H1-ii) means that there is no prey consumption as the prey level is zero and that if some prey are present, the predator are able to find and consume them. Moreover, g(x) is supposed to be increasing at the origin. As a consequence g(x) may be a Holling I or II function. g(x) might even be non-monotonic (in case of e.g. prey group defence abilities) like in Holling IV, but in this case from (H1-iii) f (x) must be of type IV too or become negative for high x values (like for logistic growth). Notice however that Holling III like functional response (i.e. those with null derivative at x = 0) fall beyond the scope of this study. No hypotheses are made on the numerical response h(x) except (H1-iv) which is verified by most classical numerical responses (but does not account for intra-predatory interference nor Allee effects for the predators). Indeed, as it will be clearer in the following, our argumentation is mainly based on the fact that g(x) is positive for positive x. Note that according to the different types of biocontrol strategies recalled in the introduction, h(x) = 0 for all x would correspond to a purely inundative biocontrol, while h(x) > 0 for positive x would take into account inoculative effects too.

We now model the augmentative biological control procedure as the T -periodic releases of biological control agents (with T > 0). Releases are, by their very nature, discrete phenomena: every time nT (n ∈ N) some predators are instantly added into the system. Let us suppose that there is a fixed rate of predators µ (i.e. number of predators per unit time) to be released in the system. Notice that such a quantity µ is a direct measure of the costs of the biological control method over some time period. Henceforth it will be referred to as well as the release rate or the biological control budget.

This results in the following: at each time period T , µT predators are added to the predator population y, yielding the discrete process ∀n ∈ N, y(nT

+ ) = y(nT ) + µT, (2)
with nT + denoting the instant right after t = nT .

Combining (1) and (2) yields the following impulsive biological control model

         ẋ = f (x) -g(x)y, ẏ = h(x)y -my, ∀n ∈ N, y(nT + ) = y(nT ) + µT. (3)
Remark 2: It is quite easy to check from (H1) that the proposed model is, as required for biological models in general, a non-negative system, i.e. the non-negative orthant of the state space is positively invariant for system (3).

Then for every non-negative initial conditions x 0 and y 0 at initial time t 0 (that only make sense with respect to the considered problem) model (3) produces non-negative trajectories x(t) and y(t) for all times t ≥ t 0 .

Special cases of model (3) with explicitly formulated functions f (x), g(x) and h(x) have been previously studied in [START_REF] Liu | The dynamical behaviors of a Lotka-Volterra predator-prey model concerning integrated pest management[END_REF][START_REF] Liu | Complex dynamics of holling type ii lotka-volterra predator-prey system with impulsive perturbations on the predator[END_REF]. Their results will be discussed after the statement and proof of our global convergence Theorem.

Global convergence

In this section we show that provided the predator release rate µ is greater than some value, augmentative biological control is able to drive any pest population to zero. This result is summarised into the following Theorem.

Theorem 1. Under Hypotheses 1, model (3) possesses a "pest free" T -periodic solution

(x p (t), y p (t)) = 0, µTe -m(t mod T ) 1 -e -mT (4)
which is locally asymptotically stable if and only if

µ > m f (0) g (0) , (5)
and globally asymptotically stable if

µ > S sup x≥0 m f (x) g(x) . ( 6
)
Remark 3: Note that from (H1-iii) S as defined in (6) does exist.

Proof: We first focus on the "pest free" set: {(x, y) ∈ R 2 + , x = 0}, which is clearly invariant by system (3) through (H1-i) and (H1-ii). Within this set, according to (H1-iv), system (3) becomes

         ẋ = 0 ẏ = -my ∀n ∈ N, y(nT + ) = y(nT ) + µT (7)
that yields: y((n + 1)T + ) = y(nT + )e -mT + µT . It is clear that the sequence (y(nT + )) n∈N has a single and globally stable equilibrium y = µT /(1e -mT ). Then, the trajectories of system [START_REF] Onofrio | Metamodeling tumorimmune system interaction, tumor evasion and immunotherapy[END_REF] globally converge towards the T -periodic solution [START_REF] Onofrio | Tumour eradication by antiangiogenic therapy: analysis and extensions of the model by hahnfeldt et al[END_REF].

We now concentrate on the stability of the pest free solution (4) for system [START_REF] Onofrio | Stability properties of pulse vaccination strategy in SEIR epidemic model[END_REF], that is to say we do not restrict ourselves to the pest free set. We change system (3) variables to consider the deviations from the pest free solution that are denoted ( x, ỹ) so that:

( x(t), ỹ(t)) = (x(t), y(t)) -(x p (t), y p (t)) that yields:    ẋ = f ( x) -g( x)( ỹ + y p (t)) ẏ = h( x)(y p (t) + ỹ) -m ỹ (8)
We first investigate local asymptotic stability (LAS) of the periodic solution (x p (t), y p (t)). Assuming ( x, ỹ) are small, we get from (8) at first order in x and ỹ

   ẋ = ( f (0) -g (0)y p (t)) x ẏ = h (0)y p (t) x -m ỹ (9)
From Theorem 1 in [START_REF] Lakmeche | Bifurcation of non trivial periodic solutions of impulsive differential equations arising chemotherapeutic treatment[END_REF], (x p (t), y p (t)) is LAS if and only if

e -T 0 mdτ < 1 and e T 0 ( f (0)-g (0)y p (τ))dτ < 1.
The former condition is trivial and the latter gives

f (0)T < T 0 g (0) µTe -mτ 1 -e -mT dτ ⇔ µ > m f (0) g (0) . ( 10
)
so that the LAS of (x p (t), y p (t)) is proved.

We now focus on the global asymptotic stability (GAS) of the pest free solution. From now on, we assume that system (3) is initiated at (x 0 , y 0 ) at time t 0 ≥ 0, i.e. system ( 8) is initiated at ( x0 , ỹ0 ) = (x 0 , y 0y p (t 0 )) at time t 0 ≥ 0.

Let us consider the function

G( x) = m x κ 1 g(s)
ds for some κ > 0. Since from (H1-ii) g(x) is positive for positive x, G(.) is an increasing function from x = 0, where it goes to -∞ since g(.) is locally Lipschitz on R + . In the following we investigate G( x(t)) behaviour as time t goes to infinity, we have

G( x(t)) -G( x0 ) = m x(t) x0 1 g(s)
ds,

= m t t 0 ẋ(τ) g( x(τ)) dτ, = t t 0 m f ( x(τ)) g( x(τ)) -m( ỹ(τ) + y p (τ)) dτ.
Considering system (8) together with (H1-iv) and Remark 2, we get ẏ ≥ -m ỹ, so that, through standard arguments ỹ(t) ≥ min(0, ỹ0 )e -m(t-t 0 ) .

Using [START_REF] Erdlenbruch | Cyclical versus non-cyclical harvesting policies in renewable resource economics[END_REF] and the definition of S in (6), we then have for all t ≥ t 0

G( x(t)) -G( x0 ) ≤ t t 0 S -my p (τ) -m min(0, ỹ0 )e -m(τ-t 0 ) dτ, = t 0 T +1 T t 0 (S -my p (τ)) dτ + t T - t 0 T -1 T 0 (S -my p (τ)) dτ + t t T T (S -my p (τ)) dτ + min(0, ỹ0 ) e -m(t-t 0 ) -1 , since y p (t) is T -periodic.
It is clear that the first, third and fourth term of the right hand side of ( 12) are upper bounded. Now, suppose that (6) holds, then in a similar way to [START_REF] Eilenberg | Suggestions for unifying the terminology in biological control[END_REF] T 0 (Smy p (τ)) dτ < 0, so that, since t T goes to infinity as t does, the right hand side of (12) goes to -∞. Then, condition [START_REF] Onofrio | Tumorimmune system interaction: Modeling the tumor-stimulated proliferation of effectors and immunotherapy[END_REF] implies

lim t→+∞ G( x(t)) = -∞ ⇔ lim t→+∞ x(t) = 0.
One can now prove that ỹ converges to zero as well. Indeed, from [START_REF] Erdlenbruch | Cyclical versus non-cyclical harvesting policies in renewable resource economics[END_REF] ỹ0 ≤ 0 ⇒ ỹ(t) ≥ ỹ0 e -m(t-t 0 ) , so that, either ỹ(t) converges to zero in infinite time from below, or it reaches the positively invariant region where ỹ ≥ 0 in finite time. Up to initial time t 0 translation, we then have to consider the positive ỹ0 only. Since x converges to zero and h(0) = 0, it is clear that there exists a time t m such that

∀t > t m , h( x(t)) ≤ m 2 ⇒ ∀t > t m , ẏ ≤ h( x)y p (t) - m 2 ỹ.
Since h( x)y p (t) goes to zero as t goes to infinity, so does ỹ. We have shown that if condition (6) holds, (0, 0) is globally attractive for (8), i.e. (x p (t), y p (t)) is globally attractive for (3).

It is then easily checked through l'Hospital rule, (H1-i) and (H1-ii) that condition [START_REF] Onofrio | Tumorimmune system interaction: Modeling the tumor-stimulated proliferation of effectors and immunotherapy[END_REF] implies that the necessary and sufficient condition (5) for LAS of (x p , y p ) holds true, so that GAS follows. 2

Some comments should be given on Theorem 1. In general, the local stability condition [START_REF] Onofrio | A general framework for modeling tumor-immune system competition and immunotherapy: Mathematical analysis and biomedical inferences[END_REF] and the global stability condition (6) are different. However there are some cases where these two conditions are similar. Indeed consider a classical Lotka Volterra (LV) system where f (x) = αx and g(x) = β x, then condition (6) becomes

µ > sup x≥0 m f (x) g(x) = mα β ,
which is also condition ( 5) applied to such a system. Then for a LV system (or systems whose prey equation is similar to the LV prey equation, see e.g. [START_REF] Tang | Models for integrated pest control and their biological implications[END_REF]) a necessary and sufficient condition for GAS of the pest-free solution ( 4) is µ > mα β . This result was previously proved for the LV model in [START_REF] Liu | The dynamical behaviors of a Lotka-Volterra predator-prey model concerning integrated pest management[END_REF]. Now consider a Rosenzweig-Mac Arthur like model with logistic growth of the prey and a Holling II functional response of the predator

f (x) = rx 1 - x K and g(x) = ax c + x . ( 12
)
Then condition (6) becomes

µ >    mr(K+c) 2 4Ka for 0 ≤ c ≤ K mrc a for c ≥ K
Since for such a system, condition (5) becomes µ > mrc a , for a carrying capacity K smaller than the half saturation constant c of the Holling II, the local stability condition is necessary and sufficient to ensure GAS of the pest free solution [START_REF] Nundloll | Two models of interfering predators in impulsive biological control[END_REF]. However if K > c, the GAS condition of the pest free solution is then stronger than the LAS condition.

If the LAS condition holds true, but the GAS does not, Theorem 1 does not allow to conclude on the global behaviour of the system. Indeed, a model verifying [START_REF] Fenton | Optimal application strategies for entomopathogenic nematodes: integrating theoretical and empirical approaches[END_REF] has been previously studied in [START_REF] Liu | Complex dynamics of holling type ii lotka-volterra predator-prey system with impulsive perturbations on the predator[END_REF] and was shown to possibly exhibit co-existing attractors: this is a possibility here and depending on the initial conditions, trajectories of the system may converge either to the pest free set or to another attractor. This conclusion obviously applies to systems for which condition ( 5) and ( 6) are different. Special cases of model (3) were also studied by [START_REF] Liu | Complex dynamics of holling type ii lotka-volterra predator-prey system with impulsive perturbations on the predator[END_REF][START_REF] Liu | The dynamical behaviors of a Lotka-Volterra predator-prey model concerning integrated pest management[END_REF] as the LAS condition of the pest free set (5) was not satisfied. The classical methodology is to prove the existence of a positive periodic solution through a bifurcation Theorem proposed by [START_REF] Lakmeche | Bifurcation of non trivial periodic solutions of impulsive differential equations arising chemotherapeutic treatment[END_REF] and then to show the existence of various routes to chaos through numerical simulations (see also [START_REF] Liu | The dynamics of a prey-dependent consumption model concerning impulsive control strategy[END_REF][START_REF] Liu | Dynamic complexities of a Holling I predator-prey model concerning periodic biological and chemical control[END_REF][START_REF] Negi | Dynamics in a Beddington-DeAngelis prey-predator system with impulsive harvesting[END_REF][START_REF] Tang | Multiple attractors in stage-structured population models with birth pulses[END_REF] for similar approaches on predator-prey systems).

In this paper however, with respect to our application we concentrate on the pest eradication problem only. In the case where condition (5) holds true, but condition (6) does not, the pest free solution is still locally stable but not globally. Such a release rate has the advantage of being smaller than the one required for GAS so that it would be cheaper for the grower to choose it as such: it allows good control of small pests invasions. The price to pay is however the risk of considerable crop damage if the system escapes to an alternative attractor to the pest free set.

It is to be noted that the stability conditions ( 5) and ( 6) on the release rate µ are both independent of the release period T . Therefore, once a release rate µ fulfilling a stability condition is chosen, pest eradication could be achieved independently of the choice of the release period T . Indeed, both infrequent and large releases (T and µT large) or frequent and small releases (T and µT small) would ultimately drive the pest population to zero. Though there is no difference in the long term, different values of the release period T might however result in different outcomes for the transient dynamics of the system. Hence in the following, we will take advantage of the adjustable parameter T to address the practically important question: how to "best" (in a sense that will be detailed in what follows) deploy a given biological control budget µ to eradicate pest outbreaks ?

Optimal Choice of the Release Period in the Preemptive Case

We first show two rather general results on one-dimensional non autonomous differential equations that are piecewise monotonous. These results will be central to the analysis we present afterwards; however the reader interested by the application only can go directly to section 3.2.

Preliminary results

Lemma 1. Consider a one-dimensional non-autonomous differential equation initiated at x(t 0 ) = x 0 ẋ = w(t),

with w(t) a T -periodic, negative, increasing function which is continuous on (0, T ) and such that

T 0 w(τ)dτ = -ωT. ( 14
)
Let x f < x 0 be the fixed terminal value of x reached at t = t f (t 0 ) and k =

x 0 -x f ωT , then either:

(i) x 0 -x f ωT
= k and (t f (t 0 )t 0 ) = kT for all t 0 .

(ii)

x 0 -x f ωT

> k and (t f (t 0 )t 0 ) is maximum at and only at the uniquely defined t * 0 such that t f (t * 0 ) = (k + 1)T .

Notice that from [START_REF] Funasaki | Invasion and chaos in a periodically pulsed mass-action chemostat[END_REF] and since w(t) is increasing in t, w(.) has necessarily a discontinuity with a negative jump at t = T .

Proof:

Without loss of generality we assume that t 0 ∈ [0, T ). Define r(t 0 ) such that t f (t 0 ) = t 0 + kT + r(t 0 ). We have from ( 13)

x f -x 0 ωT = 1 ωT t 0 +kT t 0 w(τ)dτ + 1 ωT t 0 +kT +r(t 0 ) t 0 +kT
w(τ)dτ, so that, from ( 14) and the fact that

w(t) is T -periodic k - x 0 -x f ωT = 1 ωT t 0 +r(t 0 ) t 0 w(τ)dτ. ( 15
)
From the definition of k, the left hand side of (15) belongs to (-1, 0], so that, since w(.) is negative and verifies [START_REF] Funasaki | Invasion and chaos in a periodically pulsed mass-action chemostat[END_REF], we have

r(t 0 ) ∈ [0, T ).

Suppose that

x 0 -x f ωT = k, then the left hand side of (15) cancels and r(t 0 ) = 0, so that (t f (t 0 ) -t 0 ) = kT for all t 0 .

Else, suppose that [START_REF] Grevstad | Factors influencing the chance of population establishment: Implications for release strategies in biocontrol[END_REF] with respect to t 0 and rearranging the terms we get for

x 0 -x f ωT > k, then, differentiating
t 0 0 and t 0 + r(t 0 ) T dr dt 0 = w(t 0 ) w(t 0 + r(t 0 )) -1. (16)
Since w(t) is piecewise increasing, T -periodic and r(t 0 ) ∈ [0, T ), we have

t 0 + r(t 0 ) < T ⇒ w(t 0 + r(t 0 )) > w(t 0 ), and t 0 + r(t 0 ) > T ⇒ w(t 0 + r(t 0 )) < w(t 0 ),
which with (16) yields t 0 + r(t 0 ) < T ⇒ dr dt 0 < 0, and t 0 + r(t 0 ) > T ⇒ dr dt 0 > 0.

Then as r(t 0 ) is continuous in t 0 , it has necessarily a maximum if and only if t 0 = t * 0 such that r(t * 0 ) = T -t * 0 . Thus (t f (t 0 )t 0 ) is also maximum if and only if t 0 = t * 0 and we have t f (t * 0 ) = (k + 1)T . The uniqueness of t * 0 is easily obtained from [START_REF] Grevstad | Factors influencing the chance of population establishment: Implications for release strategies in biocontrol[END_REF].

2

Lemma 1 is illustrated on figure 1. Notice that, although we do not show it here, Lemma 1 is still valid when w(t)

is locally increasing at t = 0 + and t = T -only and non-decreasing elsewhere.

Lemma 2. Consider equation [START_REF] Fenton | Evaluating the efficacy of entomopathogenic nematodes for the biological control of crop pests: A nonequilibrium approach[END_REF] and other notations as in Lemma 1. Then

max t 0 t f (t 0 ) -t 0 ≥ x 0 -x f ω
with the equality holding if and only if T =

x 0 -x f nω for some n ∈ N * .

Proof: Suppose that T =

x 0 -x f nω for some n ∈ N * . Then

n = x 0 -x f ωT and from case (i) of Lemma 1 (t f (t 0 ) -t 0 ) = nT = x 0 -x f ω
for all t 0 .

Figure 1: Graphical illustration of Lemma 1 for w(t) increasing in t. Either case (i) holds and t f (t 0 )t 0 = kT for all t 0 (thick line), or case (ii) holds and t f (t 0 )t 0 is maximum for t * 0 such that r(t * 0 ) = Tt * 0 (dashed line). It is shown in Lemma 2 that r(t 0 ) > 0 if case (i) does not hold.

Now suppose that T

x 0 -x f nω for all n ∈ N * and let k and r(t 0 ) be as in Lemma 1. Then

t f (t 0 ) -t 0 = kT + r(t 0 ),
with r(t 0 ) ∈ (0, T ), since r(t 0 ) = 0 leads to the contradiction

x 0 -x f
ωT ∈ N * . Using this in equation [START_REF] Grevstad | Factors influencing the chance of population establishment: Implications for release strategies in biocontrol[END_REF], we get at

t 0 = t * 0 x 0 -x f ω = t f (t * 0 ) -t * 0 -r(t * 0 ) - 1 ω T t * 0 w(τ)dτ, since t * 0 + r(t * 0 ) = T .
Obiously, as r(t * 0 ) < T by definition, t * 0 > 0 and as w(t) is increasing on (0, T ), we have from ( 14)

1 T -t * 0 T t * 0 w(τ)dτ > -ω, so that x 0 -x f ω < t f (t * 0 ) -t * 0 -r(t * 0 ) + T -t * 0 = t f (t * 0 ) -t * 0 2
We now come back to our main purpose, the identification of the "best" biocontrol agents release policy in model (3).

Statement of the optimisation problem and optimal choice of the release period

We focus on preemptive/preventative use of biological control agents releases: we suppose that biological control agents are released in anticipation of pests outbreaks, so that natural enemies would be able to fight the pest right at the time of their invasion. Hence in the following, we suppose that preventive releases are being performed and that the system is in the invariant set { ỹ ≥ 0} as a pest population x 0 invades the crop at some unforeseen moment t 0 .

Our goal is to provide recommendations on how to deploy such a preemptive biological control strategy (i.e. how to choose T ) to minimise crop damage due to a pest outbreak. To evaluate crop damage induced by the pests, we use the the concept of "Economic Injury Level" (EIL) that has been introduced from the early bases of theoretical biological control [START_REF] Stern | The integrated control concept[END_REF]. EIL (denoted x in the following) is defined as the "lowest (positive) pest population level that will cause economic losses on the crop". x is then a constant parameter of our problem. For a given pest outbreak

x 0 > x, we consider that the damage cost J is an increasing function of the time spent by the pest population above the EIL x, i.e.

J = t f t 0 γ(τ)dτ, (17)
with t f defined such that x(t f ) = x and γ(.) > 0. It is clear that minimising J is equivalent to minimising Π = (t ft 0 ) so that we will only concentrate on the latter in the following.

We are left with two more problems. The first one is that, for given x 0 and T , the damage cost J (or Π) strongly depends on the pest outbreak instant t 0 , i.e. for a given T the same invading pest level x 0 yields different damage costs depending on the value of t 0 . Without loss of generality, from now on we will consider that t 0 ∈ [0, T ). We seek to provide the safest biological control procedure so that we will concentrate on the damage cost for the worst t 0 , ensuring any other t 0 would result in smaller damage costs. Hence in the following, we seek the time period T that minimises the damage cost for its worst t 0 i.e. we seek to minimise max

t 0 ∈[0,T ) Π.
The remaining problem is that our analysis cannot be achieved for the general model [START_REF] Onofrio | Stability properties of pulse vaccination strategy in SEIR epidemic model[END_REF]. It is however possible to carry it out for a local approximation of ẋ about x = 0 as well as for an upper-bound of ẋ. Indeed, consider the local approximation of ẋ near (x p , y p ). From (9), and considering the change of variable z 1 = m g (0) ln x x (which is obviously increasing in x), we get ż1 = S lmy p (t), where S l = m f (0) g (0) ≤ S. Now we consider the general ẋ equation of model (3) under the hypothesis that the system is in the positively invariant set { ỹ ≥ 0}. We get for all t ≥ t 0 ẋ ≤ f (x)g(x)y p (t).

Let us make the change of variable z 2 = m x x 1 g(x) dx (which is increasing in x), we get ż2 ≤ Smy p (t), whose right hand side is similar to the ż1 equation.

In the rest of the paper we will then focus on the system ż = σmy p (t), [START_REF] Jacobson | Compatibility of beauveria bassiana (balsamo) vuillemin with amblyseius cucumeris oudemans (acarina : Phytoseiidae) to control frankliniella occidentalis pergande (thysanoptera : Thripidae) on cucumber plants[END_REF] with z standing for z 1 or z 2 and σ for S l or the original S. Through equation ( 18) we will thus study in one step the exact local approximation of our system as well as an upper bound of it, yielding a locally optimal result as well as a globally sup-optimal one.

Since the biological control program must eradicate the pest population, we assume that the stability condition of the pest free solution holds true, we have

Hypothesis 2 (H2):

The biological control budget is chosen such that: µ > σ

We first show that (see Appendix A):

Proposition 1. There exists a release period T such that if T < T , the pest population is decreasing for t ≥ t 0 .

Proposition 1 is of practical interest regarding growers concerns: indeed it is not interesting to invest in preventative releases of biological control agents that allow an invading pest population to proliferate, even temporarily, after the moment of invasion. It is much more desirable to choose a procedure that will always make the pest population decreasing after the invasion. Thus in what follows we assume that Hypothesis 3 (H3): The release period is chosen such that T ∈ (0, T ).

We now seek the moment of invasion t 0 that maximises the damage cost J (i.e. that maximises the damage time Π) for a given invading pest population x 0 > x (corresponding to z 0 > 0) and a given release period T ∈ (0, T ). Notice that with the z formalism, we have Π = t ft 0 with z(t 0 ) = z 0 and z(t f ) = 0. Since from (H2) and (H3) equation [START_REF] Jacobson | Compatibility of beauveria bassiana (balsamo) vuillemin with amblyseius cucumeris oudemans (acarina : Phytoseiidae) to control frankliniella occidentalis pergande (thysanoptera : Thripidae) on cucumber plants[END_REF] falls within the scope of Lemma 1 with ω = µσ . We easily have:

Lemma 3. Consider equation [START_REF] Jacobson | Compatibility of beauveria bassiana (balsamo) vuillemin with amblyseius cucumeris oudemans (acarina : Phytoseiidae) to control frankliniella occidentalis pergande (thysanoptera : Thripidae) on cucumber plants[END_REF], assume that z 0 > 0 and T ∈ (0, T ) are fixed, and let k = z 0 (µ-σ )T , then either:

(i) z 0 = (µσ )kT and Π = kT for all t 0 (ii) z 0 > (µσ )kT and max t 0 Π(t 0 ) = (k + 1)Tt * 0 is reached at, and only at, the uniquely defined t * 0 such that t f (t * 0 ) = (k + 1)T .

We now investigate which T ∈ (0, T ) minimises the worst case damage time max t 0 Π (i.e. that minimises the worst case damage max t 0 J). Applying Lemma 2, we get the following result:

Theorem 2. Suppose z 0 > 0 is fixed. Let T 1 = z 0 µ -σ . ( 19
)
Then there exists n 0 ∈ N * such that min T ∈(0, T ) max t 0 Π(t 0 ) = T 1 is reached at, and only at, T = T n = T 1 n for all integer n > n 0 .

Proof: Consider T 1 that is positive since µ is assumed to be greater than σ . It is then clear that there exists an integer n 0 such that for all integer n > n 0 , T n < T so that equation ( 18) verifies the condition of Lemma 1 for all n > n 0 .

The remaining part of the proof is straightforward using Lemma 2. 2

Before commenting on Theorem 2, we notice that for any T < T one can express max t 0 (Π) as a function of T 1 , T and t * 0 . Indeed, integrating (18) between t 0 and t f (t * 0 ) for some T < T yields

z(t f ) = 0 = z 0 + σ max t 0 Π - mµT 1 -e -mT T t * 0 e -mτ dτ + k T 0 e -mτ dτ = z 0 + σ max t 0 Π -(k + 1)µT + µT 1 -e -mt * 0 1 -e -mT (20)
so that, recalling from Lemma 3 that (k + 1)T = t * 0 + max t 0 Π, we get max

t 0 Π = T 1 + µ µ -σ 1 -e -mt * 0 1 -e -mT T -t 0 * (21)
Theorem 2 shows that for a given z 0 we have a countable infinite number of release periods, the T n , that solves our optimisation problem.

However, we have two difficulties with Theorem 2. On the one hand, through equation ( 18) we studied both the local optimisation problem and the global suboptimal one. To solve these two, we need to choose a release period T = T n , which depends on σ . In the local optimisation case, σ equals S l = m f (0) g (0) while in the global suboptimal one, σ equals S = sup x≥0 m f (x)

g(x)

. There is no reason why S S l would be rational, so that Theorem 2 states that it is not possible to choose a release period T that solves both the local and global optimisation problem; the best that can be done is to choose T to be as close as possible to the solution of both optimisation problems.

On the other hand, it is to be noticed that T 1 , which determines the optimal values of the release period, depends on some "biological" values, namely the invading pest population level x 0 (i.e. z 0 ) and the combination of parameters σ . This is quite a big drawback since biological parameters are usually badly known and the initial invading pest level x 0 is not known at all. Clearly such model-based optimisation approaches should take these uncertainties into account [START_REF] Loehle | Control theory and the management of ecosystems[END_REF].

Hence, in the following, we will show how T should be chosen in order to ensure that the worst case damage time max t 0 Π remains close to its minimal value, independently of the value of x 0 and σ . In doing so, we will provide an almost optimal choice of T that solves both the local and global optimisation problem and that is robust to the value of x 0 and to uncertainty of the parameters.

Robustness Analysis

In the preceding sections, we wanted our choice of T to minimise the time that it takes to reach x from x 0 for the worst possible initial time t 0 . In addition to that, we now also want this to be robust with respect to the uncertainty on the initial condition x 0 and to the model parameters that we will denote p in the rest of the paper (p is simply the pair (σ , m)). Indeed, in the poorly measured conditions of crop culture and with the accompanying scarcely known models, it is fundamental for a control method to be robust with respect to these uncertainties. We then suppose that the "true" x 0 and p belongs to some compact sets: Hypothesis 4 (H4): The uncertainties on the initial condition and the model parameters are bounded, with z 0 ∈ [z 0 , z 0 ] denoted Z and p that belongs to a compact set P.

We still assume that the invading population is greater than the EIL x (i.e. z 0 > 0) and that the stability condition holds, i.e. µ > max p∈P (σ ). Our optimisation procedure will be robust provided our choice of T minimises max t 0 Π for the worst z 0 and p choice within the sets defined in Hypotheses (H4).

Notice that the minimum (with respect to T ) of max t 0 Π is equal to T 1 that does depend on both z 0 and p. Thus in our robustness analysis it makes more sense to consider the deviation of max t 0 Π from T 1 , rather than the absolute value of max t 0 Π. Hence, we first focus on the evolution of max

(z 0 ,p) ∈ Z×P max t 0 Π(t 0 , z 0 , p) -T 1 (z 0 , p)
with respect to T . We have the following result:

Theorem 3. The function T → max (z 0 ,p) ∈ Z×P max t 0 Π(t 0 , z 0 , p) -T 1 (z 0 , p) (22)
is increasing for T smaller than some positive constant T L . Moreover

lim T → 0 + max (z 0 ,p) ∈ Z×P max t 0 Π(t 0 , z 0 , p) -T 1 (z 0 , p) = 0. (23)
Proof: From equation ( 21), we have max

(z 0 ,p) ∈ Z×P max t 0 Π(t 0 , z 0 , p) -T 1 (z 0 , p) = max p ∈ P µ µ -σ max z 0 ∈ Z 1 -e -mt * 0 (T,z 0 ,p) 1 -e -mT T -t * 0 (T, z 0 , p) .
We first concentrate on max

z 0 ∈ Z 1 -e -mt * 0 (T,z 0 ,p) 1 -e -mT T -t * 0 (T, z 0 , p) (24)
for fixed T and p. Let us introduce

T l (z 0 , z 0 , p) = min T (p), z 0 -z 0 2(µ -σ ) ,
with T (p) defined in Proposition 1.

In the following, we consider that T belongs to 0, T l (z 0 , z 0 , p) . On the one hand, this ensures through Proposition 1 that, independently of the actual parameters, the pest population (i.e. z in our new variables) always decreases after the invasion time t 0 . On the other hand, we have

z 0 -z 0 T (µ -σ ) > 2,
so that there exists n ∈ N and z 01 , z 02 ∈ [z 0 , z 0 ] such that

z 01 T (µ -σ ) = n and z 02 T (µ -σ ) = n + 1.
From Theorem 2, T is thus an optimal period for both the initial conditions z 01 and z 02 .

For all z 0 ∈ (z 01 , z 02 ), case (ii) of Lemma 3 holds for k = n , since k = z 0 T (µ-σ ) . Then, since T is an optimal period for the initial condition z 01 , we should have by continuity lim

z 0 →z + 01 Π(t * 0 (T, z 0 , p)) = (n + 1)T -lim z 0 →z + 01 t * 0 (T, z 0 , p) = T 1 (z 01 , p), so that lim z 0 →z + 01 t * 0 (T, z 0 , p) = (n + 1)T -T 1 (z 01 , p), = (n + 1)T - z 01 µ -σ T = T.
A similar argumentation on z 02 yields lim z 0 →z - 02 t * 0 (T, z 0 , p) = 0.

From [START_REF] Lakmeche | Bifurcation of non trivial periodic solutions of impulsive differential equations arising chemotherapeutic treatment[END_REF] and since max t 0 Π = (k + 1)Tt * 0 , we can conclude that t * 0 (T, z 0 , p) is a continuous function of z 0 . Then, as z 0 covers [z 01 , z 02 ], t * 0 (T, z 0 , p) reaches its whole interval of definition [0, T ] so that we have:

max z 0 ∈ [z 01 ,z 02 ]
1e -mt * 0 (T,z 0 ,p)

1 -e -mT T -t * 0 (T, z 0 , p) = max t * 0 ∈[0,T ] 1 -e -mt * 0 1 -e -mT T -t * 0 , = max z 0 ∈ Z 1 -e -mt * 0 (T,z 0 ,p)
1e -mT Tt * 0 (T, z 0 , p) .

Differentiating 1-e -mt * 0

1-e -mT Tt * 0 with respect to t * 0 , we show that it reaches its maximum for

t * 0 = 1 m ln mT 1 -e -mT . (25)
t * 0 > 0 since mT > 1e -mT . In addition, t * 0 < T since e mT > mT 1-e -mT , then max

z 0 ∈ Z 1 -e -mt * 0 (T,z 0 ,p) 1 -e -mT T -t * 0 (T, z 0 , p) = 1 -e -m t * 0 1 -e -mT T -t * 0 , = 1 m mT 1 -e -mT -1 -ln mT 1 -e -mT , H(T, p). ∂ H(T, p) ∂ T = e -mT (e mT -1 -mT ) (1 -e -mT ) 2 1 - 1 -e -mT mT > 0.
Moreover, through l'Hospital rule, we have for all p ∈ P, lim T → 0 + H(T, p) = 0.

Let us introduce:

T L = min p∈P T l (z 0 , z 0 , p) > 0.
Through the above analysis, it has been shown that for T ∈ (0, T L ) and p ∈ P, H(T, p) is an increasing function of T with lim T →0 + H(T, p) = 0.

We now come back to the function defined in [START_REF] Liu | Dynamic complexities of a Holling I predator-prey model concerning periodic biological and chemical control[END_REF]. We have max

(z 0 ,p) ∈ Z×P max t 0 Π(t 0 , z 0 , p) -T 1 (z 0 , p) = max p∈P µ µ -σ H(T, p) .
Then, from H(T, p) properties, for T smaller than T L the worst deviation of max t 0 Π from T 1 (according to (z 0 , p) ∈ Z × P) is an increasing function of T . Moreover lim

T → 0 + max (z 0 ,p) ∈ Z×P max t 0 Π(t 0 , z 0 , p) -T 1 (z 0 , p) = 0.
Π(t 0 , z 0 , p) -T 1 (z 0 , p) = 0 ( 26
)
is achieved for T = 0 + .

With the help of Theorem 3, we know that provided T is not too large, the smaller T is chosen, the less will be the worst case (according to parameters and initial condition uncertainty) deviation of max t 0 Π from its minimal value T 1 .

It has moreover been shown that as T tends to 0, the worst case deviation of max t 0 Π from its minimal value goes to 0 as well.

We conclude from this analysis that the robust choice of T , as defined by the criterion [START_REF] Loehle | Control theory and the management of ecosystems[END_REF] consists in choosing T = 0. However, it is obviously not possible to achieve such a T in practical applications: it would consist in continuously releasing predators into the crop culture (see also the discussion in Section 3.4). As an alternative, and since the argument of the infimum in ( 26) is an increasing function of T for T ∈ (0, T L ), it is then advised to take T as small as possible, i.e. to perform as small and as frequent introductions as possible, to achieve the best possible robustness.

A last remark should be made about the need for robustness. In the preceding analysis, we have concentrated on a single model in the form [START_REF] Jacobson | Compatibility of beauveria bassiana (balsamo) vuillemin with amblyseius cucumeris oudemans (acarina : Phytoseiidae) to control frankliniella occidentalis pergande (thysanoptera : Thripidae) on cucumber plants[END_REF]. We should remember that this model covers two cases: a local linearisation of the system (which yields a certain value of S) and a system that globally upper-bounds the actual nonlinear system (which yields another value of S). At the very least, the choice of T that we made should work well in both those cases. We can then conclude that, even in the unrealistic case where the initial conditions and the parameters are very well-known, one should make a choice of T that is sufficiently robust, so that both global and local approximations are covered, i.e. one should take T as small as possible.

3.4. As small and as frequent introductions as possible ?

In the preceding, it was supposed that the overall cost of the pest management program was determined by those of the biocontrol agents to be released (µ which is fixed) and of the damage J. Minimising the pest management program cost was then equivalent to minimising J. In practice, additional expenses may also be incurred by higher frequencies of releases, in terms for instance of fixed labor costs related to the time spent doing one biocontrol agents release on the crop. In this case, assuming a fixed spending for each predator release, leads to add a positive, convex and decreasing in the release period T function to the damage J to estimate the overall pest management cost. The overall worst case cost would then be the sum of a positive, convex and decreasing function in T and of a positive, convex and increasing function in T (defined on the basis of ( 22)). It is then clear that this overall worst case cost would have a minimum with respect to the release period T : below this value, the cost would be decreasing in T and above it would be increasing.

Precisely determining the value of T for which the overall worst case cost is minimum relies however on the perfect knowledge of the function [START_REF] Liu | Dynamic complexities of a Holling I predator-prey model concerning periodic biological and chemical control[END_REF], which is in turn strongly dependent on the uncertainties on the parameters and invading pest population. Hence, determining the optimal T is not really possible. However one can argue that for the high valued and sensitive crops that we consider in this paper, the labor costs related to a release remains low with respect to the costs induced by a pest outbreak. In this situation it is most probable that the optimal T would be below any reasonable biocontrol agents release period 1 . Therefore there remains fair reasons to prefer frequent and small releases over large and infrequent ones.

Numerical Simulation

In order to go beyond the worst case analysis and to illustrate the robustness result stated in Theorem 3 we performed a Monte Carlo like numerical simulation that showed how the damage time Π varied around T 1 . An impulsive predator-prey model of the form (3) was thus considered with constant a priori chosen parameters for the continuous part of the model and µ verifying condition [START_REF] Onofrio | Tumorimmune system interaction: Modeling the tumor-stimulated proliferation of effectors and immunotherapy[END_REF]. A triplet of the remaining parameters was randomly chosen: first a release period T ∈ (0, T L ), the pest invasion moment t 0 ∈ [0, T ) and the pest invasion level x 0 > x. The value of T L is easy to compute from [START_REF] Mailleret | Optimal release policy for prophylactic biological control[END_REF] since in this case the parameter set P is reduced to a singleton and we have chosen a large Z interval. For each triplet, the time to reach x, Π, as well as its deviation from T 1 , was then computed by means of numerical integration of the model. The process was repeated until a set of 2.10 5 different simulations was obtained and the results, the deviation of Π from T 1 , were plotted against the time period T . As such, we obtained a representative set of simulations of pest invasions (varying in intensity and moment of invasion) that challenged a crop under strategies of biological control varying in dose and frequency of biocontrol agents releases, as described by equation ( 2).

The results of this numerical work is presented on Figure 2, each of the dot corresponding to one of the 2.10 5

simulations. What actually corresponds to Theorem 3 is the upper envelope of the data set only (since it is the maximum with respect to parameters of the worst case, according to t 0 , of Π). Each dot corresponds to the deviation of the damage time Π from its minimal worst case T 1 for one of the 2.10 5 randomly chosen triplet T ∈ (0, T L ), t 0 ∈ (0, T ) and x 0 > x. The larger the release period T , the higher the damage time Π may be.

As it can be noticed the smaller the choice of T , the smaller the positive deviation of Π (thus of the crop damage) from its optimal value T 1 . The price to pay is however that such a choice reduces also the negative deviation of Π from T 1 . Indeed an analysis similar to the one performed in Section 3, but focusing on: min (z 0 ,p) ∈ Z×P min t 0 Π(t 0 , z 0 , p) -T 1 (z 0 , p) would show that this is a decreasing function of T with limit 0 as T tends to 0 + . Hence, if one feels lucky, a choice of T small may not be the best one since one can obtain, fortunately, smaller Π (thus smaller crop damage) if the actual t 0 is not the worst one.

Conclusion

In this contribution we studied a rather general predator-prey model with periodic impulsive additions of predators (releases) that is relevant for biological control modelling. The impulsive releases of the biological control agents were modelled on the basis of a fixed budget to be spent per time period (release rate). As such it allowed us to compare the different releases period as the same overall amount of predators was used. In a first step, it was shown that, provided the release rate of biocontrol agents was high enough, the biological control program guarantees the eradication of the prey (pest) population, whatever the value of the release periods. Then the efficiency of different release periods (for a given release rate) were compared on the assumption that the biological control program is conducted on a preventive basis. It was shown that the time needed to eradicate an unforeseen pest invasion occurring at a worst moment had a minimum that could be achieved through a countable infinite number of release periods.

However, these release period values were strongly dependent on the model parameters and a robustness analysis of this optimisation problem showed that the smaller the release period, the less the time needed to eradicate the pest invasion, for the worst invasion moment and accounting for parameter uncertainty. Hence, as an advice to biological control practitioners, our study indicates that it is less risky to release frequent small amounts of biocontrol agents rather than massive infrequent ones. This result (frequent releases of small doses yield better control) is similar to the ones of [START_REF] Yano | A simulation study of population interaction between the greenhouse whitefly trialeurodes vaporarium westwood (homoptera: Aleyrodidae), and the parasitoid encarsia formosa gahan (hymenoptera: Aphelinidae). i. description of the model[END_REF][START_REF] Yano | A simulation study of population interaction between the greenhouse whitefly trialeurodes vaporarium westwood (homoptera: Aleyrodidae), and the parasitoid encarsia formosa gahan (hymenoptera: Aphelinidae). ii. simulation analysis of population dynamics and strategy of biological control[END_REF][START_REF] Fenton | Evaluating the efficacy of entomopathogenic nematodes for the biological control of crop pests: A nonequilibrium approach[END_REF] that were obtained from simulation models. Such a biocontrol strategy is also advocated by [START_REF] De | Biological control of thrips on ornamental crops: Interactions between the predatory mite neoseiulus cucumeris (acari : Phytoseiidae) and western flower thrips, frankliniella occidentalis (thysanoptera : Thripidae), on cyclamen[END_REF] from real life experiments and by [START_REF] Skirvin | Modelling the effects of plant species on biocontrol effectiveness in ornamental nursery crops[END_REF] on the basis of a simulation model. To our knowledge, this contribution is among the firsts to recommend such a strategy on the basis of a mathematical analysis.

In a broader context, we believe that the techniques used in this work may also be of interest for optimisation problems on impulsive models from other applied fields: optimisation of vaccination strategies in epidemiology [START_REF] Shulgin | Pulse vaccination strategy in the SIR epidemic model[END_REF][START_REF] Onofrio | Stability properties of pulse vaccination strategy in SEIR epidemic model[END_REF][START_REF] Hui | Impulsive vaccination of SIR epidemic models with nonlinear incidence rates[END_REF] or optimisation of chemotherapy in cancer treatments [START_REF] Lakmeche | Bifurcation of non trivial periodic solutions of impulsive differential equations arising chemotherapeutic treatment[END_REF][START_REF] Onofrio | Tumour eradication by antiangiogenic therapy: analysis and extensions of the model by hahnfeldt et al[END_REF][START_REF] Onofrio | A general framework for modeling tumor-immune system competition and immunotherapy: Mathematical analysis and biomedical inferences[END_REF]. For instance, a first application may be to re-investigate the work by Pugliese and Gandolfi in immunology [START_REF] Pugliese | A simple model of pathogen-immune dynamics including specific and non-specific immunity[END_REF], that shows, without the explicit use of impulsive modelling techniques, that frequent host re-infections and/or vaccinations are necessary to the persistence of host immunity.

2Corollary 1 .

 1 One can deduce the following Corollary from Theorem 3 (see Appendix B): Let T = min p ∈ P T (p), then inf T ∈ (0,T ) max (z 0 ,p) ∈ Z×P max t 0 ∈ [0,T ]

Figure 2 :

 2 Figure 2: Graphical illustration of Theorem 3. Each dot corresponds to the deviation of the damage time Π from its minimal worst case T 1 for one

Indeed, and as in[START_REF] Grevstad | Factors influencing the chance of population establishment: Implications for release strategies in biocontrol[END_REF] that reported an analysis on the optimisation of the spatial deployment of a fixed amount of biological control agents, "as small as possible" (size and period of releases) should be interpreted loosely: too small releases would actually drive the system beyond the validity domain of the spatially implicit model (3).

Acknowledgements: This work was supported by the COLOR action "LutIns&co" funded by INRIA Sophia Antipolis Méditerranée and by the Project "IA2L" funded by the SPE Department of INRA. The authors are grateful to Sapna Nundloll for polishing the English.

Consider equation [START_REF] Jacobson | Compatibility of beauveria bassiana (balsamo) vuillemin with amblyseius cucumeris oudemans (acarina : Phytoseiidae) to control frankliniella occidentalis pergande (thysanoptera : Thripidae) on cucumber plants[END_REF]. Provided min t y p (t) > σ m , z is decreasing for t ≥ t 0 ; then so does the pest population x. Notice that min t y p (t) = µT e mT -1 is a decreasing function of T with lim

since the stability condition holds.

Then Proposition 1 holds true with T solution of

2

B. Proof of Corollary 1

From Theorem 3:

and is achieved for T = 0 + since the argument of the infimum in [START_REF] Mailleret | A note on semi-discrete modelling in the life sciences[END_REF] is increasing in T and tends to 0 as T does.

To complete the proof of Corollary 1 we have to study max

for T ∈ T L , T . Notice that from T L definition, we have:

so that the set T L , T is non-empty if and only if T > min (z 0 ,p)∈Z×P z 0 -z 0 2(µ-σ ) . Suppose that this holds and that T ∈ T L , T , then we have max

since otherwise it would be required that ∀z 0 ∈ Z there would exist n ∈ N such that

1e -mt * 0 (T,z 0 ,p)

1e -mT t * 0 (T, z 0 , p) , this latter being positive, which concludes the proof. 2


Vilim Filipović
email: vfilipovic@agr.hr

Joel Fábio

Kochem Mallmann

Yves Coquet

Jirka Šimůnek

Numerical simulation of water flow in tile and mole drainage systems

Keywords: Tile drainage, Mole drainage, Numerical simulation, Heavy soil, Threedimensional modeling

Tile drainage systems are sometimes not sufficient to provide favorable unsaturated conditions in the rootzone. These drainage systems then need to be supplemented with an additional high conductivity material in the trenches above the tiles or by implementing mole 2 drainage. The HYDRUS (2D/3D) model was used to evaluate the impact of such additional measures for heavy clay soil. Three types of drainage systems were simulated: i) tile drains, ii) tile drains with gravel trenches, and iii) tile drains with gravel trenches and mole drains, using either two-dimensional (the former two systems) or three-dimensional (the latter one) transport domains. Three scenarios were considered to test the efficiency of each system: i) time to drain an initially saturated system, ii) high intensity rainfall, and iii) a real case scenario. Different horizontal spacings between tile drains with or without gravel trenches were also compared with the system which included mole drainage. The results showed that the drainage system that included mole drains and gravel trenches was the most efficient. This system provided the largest drainage rate, was the first to reach steady-state in the time to drain scenario, and also efficiently reduced surface runoff. Adding mole drains to a system with tile drains and gravel trenches resulted in a large reduction of surface runoff (75%).

Simulations showed that the spacing of tile drains with or without gravel trenches would have to be 40% or 55% smaller, respectively, in order to reproduce the same water table levels as those observed for the drainage system with mole drains. Therefore, introducing mole drains in drainage systems is an efficient practice for reducing waterlogging and runoff.

Introduction

Soil drainage systems aim at limiting saturated conditions in the soil profile that can arise due to hydrological processes from above (downward water percolation) and from below (elevated groundwater). The main goal of drainage systems is to remove excess water and maintain favorable unsaturated conditions in the rootzone. The most common agricultural drainage system consists of perforated PVC tile drains installed in soil at various spacing's and depths, depending on soil hydraulic properties, climatic conditions, and cultivated crops.

When designing drainage systems, it is crucial to increase the hydraulic functioning of the entire system. These systems are mostly installed in soils with high clay content, i.e., heavy soils with very low hydraulic permeability [START_REF] Tuli | Comparison of air and water permeability between disturbed and undisturbed soils[END_REF]. In such soils, subsurface drainage systems can substantially reduce surface runoff, shorten periods of surface ponding, and lower water table [START_REF] Konyha | Effects of drainage and water-management practices on hydrology[END_REF][START_REF] Skaggs | Hydrologic and water quality impact of agricultural drainage[END_REF]. Under certain conditions, including high groundwater table, large and intensive precipitation, and heavy-textured soils, the installation of tile drains may not be sufficient to provide favorable conditions for growing crops. Additional measures may then be needed, such as using a backfill material (gravel) with a high hydraulic conductivity above tile drains or performing mole drainage. The presence of a gravel layer above tile drains up to the tilled layer may increase the efficiency of the entire system by promoting by-pass flow of water from the tilled layer directly into the drains. Mole drains are closely-spaced unlined channels of limited duration that are formed in clay subsoil using a ripper blade with a cylindrical foot, often with an expander, which helps to compact and stabilize the channel walls. Mole drainage has been recommended for heavy soils with low permeability, which would otherwise require a small drain spacing [START_REF] Hudson | The Draining of Farm Lands[END_REF]. During the construction of the mole channels the foot and expander create a 'leg slot' directly above the mole as well as fissures in the soil upwards from the mole towards the plough layer. These cracks and fissures can promote preferential flow towards mole drains [START_REF] Leeds-Harrison | Water Flow to Mole Drains[END_REF]. Mole drains are intended to improve lateral flow to tile drains and are usually used only in soils with a high clay content (>35%). If established properly and under the right conditions, they may still be operating after five years [START_REF] Harris | Effect of mole submergence on the life of mole channels[END_REF]. Tile and mole drains are nonconductive under unsaturated conditions because positive pressure must occur before water can start flowing into them [START_REF] Stormont | Impact of Unsaturated Flow on Pavement Edgedrain Performance[END_REF]. In recent years, numerical models have been developed to simulate water flow in tile-drained soils. One-dimensional (1D) models (e.g. MACRO, [START_REF] Larsbo | MACRO 5.0. A model of water flow and solute transport in macroporous soil[END_REF] use an approximate analytical solution of Darcy's equation to account for water loss through the drainage system (e.g., Hooghout, Ernst, or Boussinesq's equations), while two-or threedimensional (2D or 3D) models (e.g. HYDRUS (2D/3D), Šimůnek et al., 2006) . An example of such model is DrainMod, which has been used in many applications. For example, [START_REF] Skaggs | DrainMod: Model use, calibration and validation[END_REF] used DrainMod to simulate water flow in a subsurface-drained agricultural field in eastern North Carolina. The performance statistics indicated that the model with calibrated input data accurately predicted daily water table depths, daily drainage rates, and monthly drainage volumes. [START_REF] Singh | Calibration and validation of DRAINMOD to design subsurface drainage systems for Iowa's tile landscapes[END_REF] calibrated and tested DrainMod on two types of soils in Iowa and used it to simulate the impacts of different designs of subsurface drainage systems.

Simulation results suggested that a drainage system designed for a drainage intensity of 4.6 mm d -1 with a drain depth of 1.05 m and a drain spacing of 25 m was sufficient to maximize crop production under the prevailing local agricultural conditions.

Water flow and solute transport to tile drains has been also evaluated using the HYDRUS family of codes (Šimůnek et al., 2008). For example, the HYDRUS-2D software package was used by De [START_REF] De Vos | Nitrate leaching in a tile-drained silt loam soil[END_REF][START_REF] De Vos | Chloride transport in a recently reclaimed Dutch polder[END_REF] to simulate nitrate transport in a tile-drained layered silt loam soil in a reclaimed Dutch polder, or by [START_REF] Castanheira | Impact of mole drains on salinity of a vertisoil under irrigation[END_REF] to evaluate the impact of mole drains on salinity of a vertisol under irrigation.

Numerical modeling of drainage systems involving both tile and mole drains is much more limited. [START_REF] Snow | Predicting soil water, tile drainage, and runoff in a mole-tile drained soil[END_REF] used the APSIM-SWIM model to predict drainage rates and runoff in a mole-tile drained silty loam soil of New Zealand. There was an excellent agreement between simulated and measured drainage, as well as a reasonable agreement between measured and simulated cumulative surface runoff. Armstrong et al. (2000) performed a modeling study on a macroporous clay soil with mole and tile drains, in which they compared four "preferential flow" models (MACRO, CRACK-NP, SIMULAT and PLM) in their ability to simulate isoproturon leaching. MACRO model gave the best results, although globally the simulations showed the difficulty of deriving adequate parameters, even where relatively complete soil physical data were available. A similar study was performed by [START_REF] Besien | Simulation of water movement and isoproturon behaviour in a heavy clay soil using the MACRO model[END_REF] on a structured heavy clay soil with tile (0.75 m depth and 50 m spacing) and mole drains (0.5 m depth and 3 m spacing) in which the MACRO model was calibrated and used to investigate the leaching of isoproturon. [START_REF] Madvar | Modeling hydraulic performance of mole drains and combined mole-pipe drains[END_REF] used the SEEP/W model to simulate the hydraulic performance of mole and tile drains in heavy-textured soil. Based on their numerical study the combination of mole and tile drains resulted in economical and hydraulic improvement of the agricultural system.

Most of modeling studies considered water flow in one-or two-dimensional soil profiles, either taking into account only tile drainage without mole drainage or mole drainage without tile drainage (e.g., [START_REF] Castanheira | Impact of mole drains on salinity of a vertisoil under irrigation[END_REF], but not both in a fully threedimensional system. This is mainly because mole drains, being installed in a perpendicular direction to tile drains, make the entire system three-dimensional. Since we could not identify any study that simulated water flow in systems with tile and mole drains as a threedimensional problem, we conducted 3D numerical experiments to evaluate the performance of these systems. In our numerical simulations we considered environmental conditions that are typical for the eastern part of Croatia. Hydromorphic soils in Croatia cover an area of 1,618,500 ha, which is approximately one third (29%) of the total area occupied by agricultural soils [START_REF] Husnjak | The Croatian agricultural soils and the requirements for the amelioration measures[END_REF]. Subsurface tile drainage systems are installed on 161,530 ha [START_REF] Petošić | Drainage outflow in different pipedrainage variants on gleyic podzoluvisol in the Sava river valleyy[END_REF], with mole drainage often used to improve drainage efficiency.

The main objective of this paper was to numerically evaluate the performance of different subsurface drainage systems of increasing complexity in hydromorphic soils using a threedimensional model that can explicitly account for any drainage system. (i) The simplest system to be considered consists of a layered soil profile with parallel tile drains. Such system can be simulated using a two-dimensional simulation domain that is perpendicular to tile drains. (ii) In the second set of simulations it is assumed that a backfill material (gravel) of higher hydraulic conductivity than the original soil is placed above the drains. Such system can also be evaluated assuming a two-dimensional simulation domain perpendicular to drains.

(iii) Finally, the most complex system consists of tile drains with backfilled gravel above drains, combined with perpendicular mole drains. Such system has to be analyzed using a fully three-dimensional flow model.

Materials and methods

Theory

Subsurface water flow in the unsaturated and saturated zones is governed by the mass balance equation and the Darcy-Buckingham's law:

(1)

(2) respectively, where θ is the volumetric water

content [L 3 L -3 ], J wi is Darcy's flux [LT -1 ] in the i-th direction, S is a sink/source term [T -1 ], K(h) is the saturated/unsaturated soil hydraulic conductivity function [LT -1 ], K ij
A are the components of the dimensionless anisotropy tensor for hydraulic conductivity K A [-], H is the total hydraulic head [L] defined as the sum of the pressure head and the gravitational head, H = h+z, x i are the spatial coordinates [L], with x and y being horizontal coordinates and z the vertical coordinate directed upwards, and t is time [T]. The governing flow equation for these conditions, resulting from combining (1) and

(2), is given by the following modified form of the Richards' equation:

(3)

The unsaturated hydraulic conductivity function K(h) is given by:

K(h, x, y, z)=K s (x, y, z)K r (h, x, y, z) (4)
where K r is the relative hydraulic conductivity [-] and K s is the saturated hydraulic conductivity [LT -1 ].

Soil hydraulic functions were described using the van Genuchten -Mualem model (van Genuchten, 1980), which is defined as follows:

for 0 (5) for 0

(6) (7) 1 (8)
where θ r and θ s denote residual and saturated volumetric water contents [L 3 L -3 ], respectively;

S e is the effective saturation [-], α [L -1 ], and n [-] are retention curve shape factors, and l is a pore connectivity parameter [-]. A numerical model for simulating water flow and solute transport in two-and three-dimensional transport domains, HYDRUS (2D/3D) (Šimůnek et al., 2006), was used it this study.

Numerical Experiments and Input Data

Soil Hydraulic Properties

Numerical experiments were set up for a Eutric Haplic Calcaric Siltic Gleysol (horizons: Ap; Bg; Cr; Cg), which is located in eastern Croatia (45°09' N and 18°42' E). Disturbed soil samples were taken from each horizon and particle size distribution was determined using the pipette method [START_REF] Gee | Particle size analysis[END_REF]. Undisturbed soil samples (100 cm 3 ) were used to measure bulk density and soil hydraulic properties. The saturated hydraulic conductivity, K s , was measured using the constant head method [START_REF] Klute | Hydraulic conductivity and diffusivity: laboratory methods[END_REF]. The saturated water content, θ s , was measured using the ISO 11274:1998 sandbox method [START_REF] Clement | A simple and reliable tension table[END_REF]. Data points of the soil water retention curve were measured using the pressure plate apparatus [START_REF] Dane | Water retention and storage: Laboratory[END_REF] for applied pressures of 3, 10, 33, 100, 625, and 1500 kPa. While the saturated water content, θ s , was measured, the remaining parameters of the soil water retention curve (5) (θ r , α, and n) were optimized using the RETC software [START_REF] Van Genuchten | The RETC Code for Quantifying the Hydraulic Functions of Unsaturated Soils[END_REF] by fitting measured data. R 2 values for soil water retention curve fitting were in range from 0.91 to 0.99. The pore connectivity parameter, l, was assumed to be equal to an average value for many soils (l=0.5) [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF]. The hydraulic properties of the material that was assumed to be backfilled into the tile trenches were taken from the soil catalog of [START_REF] Carsel | Developing joint probability distributions of soil water retention characteristics[END_REF] for sand, with an increased K s value to imitate the hydraulic properties of gravel (Table 1). As only one seepage face boundary condition could be considered in HYDRUS (2D/3D), mole drains were simulated as a gravel material with hydraulic properties identical to those of the material used for tile trenches.

Drainage systems

Three different drainage systems were evaluated in our numerical simulations. The simplest system (System 1) consisted of parallel tile drains and no other additional measures (denoted below in graphs and tables as "tile"). In simulations for System 2 (denoted as "tile_gravel"), it was assumed that a backfill material (gravel) of higher hydraulic conductivity was placed above drains. Finally, the most complex system (System 3) consisted of tile drains with a backfilled gravel above drains combined with perpendicular mole drains (denoted as "tile_gravel_mole"). Fig. 1 shows a schematic of the simulated drainage systems. System 1 and 2 were evaluated assuming a two-dimensional transport domain perpendicular to tile drains (Fig. 1a), while System 3 needed a fully three-dimensional water flow model (Fig. 1b).

Transport Domain

Because of the symmetry of water flow between two parallel tile drains, it was possible to consider only a half domain between two tile drains. A two-dimensional transport domain for Systems 1 and 2 (Fig. 1a) thus had a width of 750 cm (tile drain spacing was assumed to be 15 m) and a height of 200 cm. The 8 cm diameter tile drain, which was located on the left side of the transport domain, was assumed to be centered at a depth of 96 cm from the bottom of the soil profile. The soil profile was divided into 4 soil horizons (Table 1) based on a soil survey of the area. In simulations of System 2, an additional vertical layer of gravel, representing a backfill material (material 5; Table 1), was considered. The three-dimensional transport domain for System 3 (Fig. 1b) with tile and mole drains and a backfill material above the tile drains had the same first two dimensions as the two-dimensional transport domains, i.e., a width of 750 cm and a depth of 200 cm. The third dimension, perpendicular to tile drains, was equal to 100 cm (mole drain spacing was assumed to be 2 m). The mole drain, which was perpendicular to tile drains, was located at the front of the three-dimensional transport domain (Fig. 1b) and was assumed to have a diameter of 6.5 cm, an inclination of 0.5%, and was 40 cm deep (from the soil surface) above the tile drain.

Simulated scenarios

Numerical simulations were carried out for each drainage system for three different scenarios with different initial and boundary conditions. In Scenario 1, a "time to drain" scenario was used to compare the drainage response of the three simulated systems, i.e., tile drains, tile drains with gravel trenches, and tile drains with gravel trenches and mole drains, starting from wet soil conditions. The initial pressure head was assumed to be at equilibrium conditions with a positive pressure head of 175 cm at the bottom of the transport domain which left the first 25 cm (corresponding to the tilled layer) unsaturated. No flow boundary condition was considered for all boundaries except for the small semicircle on the left of the transport domain which represented the tile drain (Fig. 1a). A seepage face boundary condition was used to represent the tile drain. Simulation time was equal to 2,400 hours (100 days).

However, as there was still some limited outflow after 2,400 hours, the time when drainage was assumed to stop was defined as the time when the drainage rate decreased below 0.001 mm h -1 . Simulations for Systems 1 and 2 were also performed in 3D to check if the results

were identical to those provided by 2D modeling.

In the "rainfall" scenario (scenario 2), the initial conditions were set up as the pressure head distribution from the end of the first scenario (steady state conditions), and a rainfall event (with an intensity of 2 mm/h and a duration set up to get a total of either 50 mm, 100 mm, or 150 mm) was applied starting at the first day of simulation. In addition to a seepage face boundary condition representing a tile drain, an atmospheric boundary condition was applied at the top of the transport domain to include rainfall. Evaporation was neglected in this scenario.

In the "real case" scenario (scenario 3), 4-year simulations (2009-2012) were performed with meteorological data collected at the Gradište station (45°09' N and 18°42' E) using the same initial conditions as in scenario 2. An atmospheric boundary condition was applied at the top of the transport domain using evapotranspiration values calculated by the Penman-Monteith approach [START_REF] Monteith | Evaporation and surface temperature[END_REF]. Transpiration was calculated formaize (Zea Mays L., 2009[START_REF] Castanheira | Impact of mole drains on salinity of a vertisoil under irrigation[END_REF][START_REF] Skaggs | DrainMod: Model use, calibration and validation[END_REF] and barley (Hordeum vulgare L., 2011) which were grown at the actual location.

Again, a seepage face boundary condition was selected to represent the tile drain.

Results and discussion

"Time to drain" scenario

To evaluate the specific drainage response of each of the three drainage systems, a "time to drain" scenario was performed. The ending time of drainage was reached at 681 h (28.4 days) for the tile system, 594 h (24.7 days) for the tile_gravel_system, and 506 h (21.1 days) for the tile_gravel_mole system. The value corresponding to 90% of the final drainage volume from each system was selected as the time when a quasi steady-state condition was reached.

According to this criterion, the systems reached quasi steady-state after cumulated outflow amounts of 4.86 mm, 6.57 mm, 9.18 mm for the tile, tile_gravel, and tile_gravel_mole systems, which corresponded to days 11, 7.1, and 6, respectively (Fig. 2). The difference in time response of the three systems may have a large influence on the crop root system, i.e., the faster the soil is drained the less consequences on the crop root system and yield can be expected [START_REF] Oosterbaan | Agricultural Land Drainage: a wider application through caution and restraint[END_REF]. The soil horizons had a relatively low permeability and, in the absence of a highly conductive layer, the drainage is very slow. When a gravel trench above the tile is added, the time response is faster (Fig. 2) since the trench is filled up with a material that has a high saturated hydraulic conductivity. The response of the system with moles is fastest since it includes both the trench and the mole drain constructed from the same highly conductive material as the trenches. Since the domain was initially saturated well above the mole drain depth, the response of the system is increasing in the order tile < tile_gravel < tile_gravel_mole.

In addition to the time response of each system, a large difference in the cumulative amount of drainage was found. This difference is due to the substitution of soil materials having different water retention properties within the simulation domain. In System 2, parts of the volumes of materials 1, 2, and 3 are replaced by material 5 to represent the gravel trench. In System 3, in addition to the trench, a part of material 2 is replaced by material 5 to represent the mole drain. The difference in water storage between System 1 and System 2 is 1.86 mm and 4.86 mm between System 1 and 3. Material 5 and materials 1, 2, and 3 have very different water retention properties, which explain the different cumulative outflows (Table 2).

Inserting a gravel trench above the tile drain increased outflow by 1.86 mm, 90% of it coming from the water initially stored in the trench itself. Adding a mole drain (System 3) further increases the outflow by 3 mm, with more than half of it coming from the gravel material representing the mole drain. From these results one can already expect a significant effect of the mole drain on the recession of the water table.

Fig. 3 shows the simulated cross sections of all three systems during the first 24 hours. The results are displayed after 2 and 24 hours. Pressure head values are lower in the tile_gravel system than in the tile system and the lowest in the tile_gravel_mole system. After 2 h of drainage, the influence of the trench in System 2 is limited to the right side of the domain while in System 3 the effect is expanding more to the left side because of the presence of the mole drain. After 24 h, there is a large drop in the water table depth for all three systems in the order (high to low) tile > tile_gravel > tile_gravel_mole. The mean pressure head after 24 hours in the domain above the tile depth (from the surface down to 100 cm) was -75 cm, -78 cm, and -84 cm for the tile, tile_gravel, and tile_gravel_mole system, respectively. These results indicate the higher efficiency of the tile_gravel_mole system in reducing the water table height compared to the tile_gravel system.

High intensity rainfall scenario

High intensity rainfall events with an intensity of 2 mm per hour and various durations, to get total amounts of 50 mm, 100 mm, or 150 mm, were applied in the second scenario. The initial condition was set equal to the final pressure head distribution at the end of the first scenario, i.e., the surface pressure head was -108 cm and the bottom pressure head was 92 cm on the right side of the domain. The cumulative outflow increased with the rainfall intensity in all three systems (Fig. 4). In this scenario, the trench and mole drain were not initially saturated as in the previous scenario, which led to the slower response time at the beginning of the experiment. The gravel material (representing the trench and mole) was not initially saturated and the time needed for its saturation induced a delay in the system response, which can be best seen for the tile_gravel_mole system during the first 4 hours (Fig. 4). There is no outflow during the first 4 hours in System 3 since the mole drain needs to be saturated to trigger drainage. The tile system without a gravel layer became almost instantaneously saturated since the rest of the domain had a very low permeability (5 mm h -1 ).

The low permeability of the soil promotes runoff during high intensity rainfalls. In Fig. 5, the total amount of surface runoff is presented for each system. One can see that the cumulative amount of runoff increases according to rainfall intensity, while the opposite was found for cumulative outflow (Fig. 4). The largest runoff occurs in the tile system, then in the tile_gravel system, and the smallest in the tile_gravel_mole system. The mole drain was very efficient in reducing surface runoff. These results are in agreement with several studies, in which authors found that artificial drainage increases the amount of infiltrating water and reduces runoff (e.g., [START_REF] Bengtson | Agricultural drainage and water quality in Mississippi delta[END_REF].

Mass balance information for all three system and three rainfall events is shown in Table 3.

Table 3 shows that the effectiveness of mole drainage in reducing surface runoff is not only due to a better efficiency in promoting drainage outflow, but also due to a better capacity for water storage than for the other two drainage systems.

Real case scenario

For the real case scenario, agroecological conditions that correspond to East Croatia were taken into account to be able to explain the behavior of each system in real field conditions. This effect is due to the small rainfall amounts at the beginning of the simulated time period which were drained by the tile_gravel system, but not by the tile_gravel_mole system, because mole drains did not become saturated enough to trigger outflow. Long term simulations showed the opposite effect on the 1 st of June 2010 when a high intensity rainfall event (72.8 mm) saturated mole drains and promoted a large outflow from the tile_gravel_mole system, while in the same time other two systems generated low outflow and large surface runoff (Fig. 6b). Adding gravel in the trench above a drain resulted in small decrease of surface runoff by (1%) However, inserting a mole drain into the system decreased the runoff by an additional 75%. The tile drain system with gravel and mole drains effectively reduced surface runoff in real field conditions. In such situation, the tile_gravel system cannot conduct such large amounts of water as the system with moles, which leads to waterlogging conditions at the surface or the presence of a high water table.

Drain spacing

Finally, simulations with different tile drain spacings were performed in order to compare the water table variations for each system and to try to obtain, either with only tile drains or tile drains with backfilled trenches, the same water table level as for System 3. The simulation set up e.g., initial and boundary conditions was taken from the "time to drain" scenario. The tile with gravel trenches and mole system had the same spacing of 15 m as in above studied scenarios. The tile and tile_gravel systems were set up with various spacings from 6 to 12 m.

Figs. 7a and 7b show the water table elevation during 96 h in the middle between two tile drains for the three systems (numbers 6, 7, 8, 9, 10, 12 corresponds to the drain spacing in meters). At the beginning, there is a rapid decrease in water table depths in all three systems, clearly visible during the first 24 hours. After that, the decrease is smaller and starts to stabilize as time approaches the end of simulation. The tile_gravel system was more efficient than the tile system at each spacing in lowering the water table (by ~2 cm). To reproduce the same water table level after 96 h (-100.2 cm) as in the tile_gravel_mole system, one would need to lower the spacing of the tile_gravel system to 9 m, or to 7 m if the system without gravel is used. System 1 and 2 with a 6-m spacing were able to lower the water table more efficiently than System 3 with a 15-m tile drain spacing and a 2-m mole drain spacing. However, a system with such small tile drain spacing is not economically sustainable. A small tile drain spacing is mostly used only in golf course construction or in civil engineering. Since the construction of mole drains requires only a simple equipment (i.e., a special kind of plough) [START_REF] Hopkins | Managing wet soils: mole drainage[END_REF], this system could significantly lower the overall cost of the drainage system.

Conclusions

Numerical experiments have been performed for three different drainage systems: i) tile drains, ii) tile drains with gravel trenches, and iii) tile drains with gravel trenches and mole drains. The tile drain system showed in all scenarios the slowest outflow response and the smallest cumulative outflow. When compared with the other two systems, the drainage system with tile drains, gravel trenches, and mole drains provides the largest efficiency.

Under real case scenario, the system with tile drains and gravel trenches had the fastest reaction time and provided drainage outflow for small intensity rainfalls. However, the system with tile drains, gravel trenches and mole drains was able to provide larger outflow during higher intensity rainfall (>15 mm). Adding mole drains to a system with tile drains and gravel trenches resulted in a very large reduction of surface runoff (75%), which can have a very important effect during high intensity rainfalls. Since the use of the drainage systems with mole drains is still limited to special agro ecological conditions, there is only limited extensive field data available which allow to test fully three-dimensional simulations of these particular drainage systems. Research needs to be expanded through field trials in order to have more information on water flow dynamics in such complex systems.

Armstrong, A., Aden, K., Amraoui, N., Diekkrüger B., Jarvis, N., Mouvet, C., Nicholls, P., Wittwer, C., 2000. Comparison of the performance of pesticide-leaching models on a cracking clay soil: results using the Brimstone Farm dataset. Agr. Water Manage. 44, 85-104.

List of Figures

Fig. 1. Transport domains with 4 or 5 soil materials considered in the numerical simulations.

(a) a two-dimensional transport domain was considered for the tile drains and the tile drains plus gravel trenches systems, (b) a three-dimensional transport domain was considered for the tile drains plus gravel trench plus mole drains system.

  use the Richards' equation to explicitly model the dynamics of the water table. In 1D models, the flow to tile drains is implemented as a sink term in the mass balance equation. The Houghoudt's equation is based on the Dupuit-Forschheimer assumptions with corrections for convergence of radial flow near the tiles. When combining the Houghoudt's equation (for the saturated zone) with the Richards equation (for the unsaturated zone), the Houghoudt's equation gives instant lateral fluxes to tile drains, while the Richards equation gives transient fluxes in the unsaturated zone, leading to a rise or drop of the water table

Fig. 6a shows

  Fig. 6a shows the cumulative outflow during 2009-2012 period. At the beginning, one can see

Fig. 2 .

 2 Fig. 2. Cumulative outflow rates for the "time to drain" scenario.

Fig. 3 .

 3 Fig. 3. Pressure head distributions after 2 and 24 hours for the tile drains, tile drains with

Fig. 4 .

 4 Fig. 4. Cumulative outflow in the second scenario with high intensity rainfall events of 50

Fig. 5 .

 5 Fig. 5. Cumulative surface runoff in the second scenario with high intensity rainfall events of

Fig. 6 .

 6 Fig. 6. Cumulative outflow (a), cumulative runoff (b) and cumulative surface flux (c) for the

Table 1 .

 1 ). van Genuchten soil hydraulic parameters for all horizons of the soil profile used in the numerical simulations.

	1

θr (cm 3 cm -3 ) θs (cm 3 cm -3 )  (cm -1 )

						n	Ks (cm day -1 )	l
	1	0-30	0.095	0.42	0.00136	1.20	12	0.5
	2	30-70	0.095	0.41	0.00212	1.18	14	0.5
	3	70-100	0.095	0.41	0.00136	1.20	12	0.5
	4	100-200	0.102	0.50	0.0121	1.41	12.8	0.5
	5	Gravel (trench and mole) material	0.005	0.42	0.1	2.1	3000	0.5

Table 2 .

 2 Water balance (ΔWSchange in storage) for the various material layers composing the three simulated systems at the end of "time to drain" scenario (mm).

			All Materials	Mat 1	Mat 2	Mat 3	Mat 4	Mat 5
	ΔWS in each system	Tile	-5.35	-1.12	-0.87	-0.68	-2.68	0.00
		Tile_gravel	-7.21	-1.14	-0.96	-0.72	-2.68	-1.71
		Tile_gravel_mole	-10.21	-1.14	-1.91	-0.72	-2.68	-3.76
	Difference between	Tile vs Tile_gravel	-1.86	-0.02	-0.09	-0.04	0.00	-1.71
	the systems	Tile vs Tile-gravel_mole	-4.86	-0.02	-1.04	-0.04	0.00	-3.76

Table 3 .

 3 Mass balance for all three simulated systems for scenario 2subscripts 50, 100, and 150 correspond to the cumulative rainfall amount (mm).

		System 150 System 250 System 350 System 1100	System 2100 System 3100 System 1150 System 2150 System 3150
	Rainfall	50.00	50.00	50.00	100.00	100.00	100.00	150.00	150.00	150.00
	Outflow	16.57	19.58	28.13	29.71	35.65	51.22	42.81	51.69	74.26
	Runoff	31.93	29.69	6.36	68.66	63.52	13.61	105.39	97.35	20.85
	ΔWS	1.49	0.72	15.51	1.63	0.83	35.18	1.79	0.96	54.89

List of Tables

 Table 1.van Genuchten soil hydraulic parameters for all materials of the soil profile used in the numerical simulations.

Table 2. Water balance (ΔWSchange in storage) for the various material layers composing the three simulated systems at the end of "time to drain" scenario (mm). Table 3. Mass balance for all three simulated systems for scenario 2subscripts 50, 100, and 150 correspond to the cumulative rainfall amount (mm).