Michael Herbst
email: m.herbst@fz-juelich.de

W Fialkiewicz
email: wieslavfialkiewicz@brgm.fr

T Chen
email: taochen@hmg.inpg.fr

Thomas Pütz

D Thiery
email: d.thiery@brgm.fr

Christophe Mouvet
email: c.mouvet@brgm.fr

G Vachaud

Harry Vereecken
email: h.vereecken@fz-juelich.de

Unsaturated methabenzthiazuron transport in a cropped lysimeter: comparison between models and measurements

The vertical water flow, heat flow and transport of the herbicide methabenzthiazuron were monitored for 627 days in lysimeters sampled at a field site close to the research centre Jülich, Germany. During this period the lysimeters were cropped with winter wheat, winter barley and oat. The models TRACE, MARTHE, ANSWERS and MACRO were applied to the lysimeter data with the scope of upscaling local scale process understanding for regional scale. MARTHE and TRACE solve the 3-d Richards' equation for variably saturated water flow. MACRO is a 1-d model based on the Richards' Equation and accounting for preferential flow in the unsaturated zone, while ANSWERS is a regional scale capacity based watershed model. Measurements of soil moisture, evapotranspiration, drainage, soil temperature, pesticide residues and leaching are used for comparison with model results.

Although the adopted models differ in terms of model concepts, the use of model performance indices proved a proper simulation of water flow for all models. The heat flow is also well described with ANSWERS, MARTHE and MACRO. Larger deviations were found between model results and measured pesticide transport. An inadequate reproduction of the measured MBT degradation was found for the available model input parameters. A very small amount of MBT leaching, observed in the measurements, was only reproduced with MACRO after strong calibration. In other respects only plant parameters were calibrated. Calibration of the crop conversion factor used for scaling of the potential evapotranspiration was found to be a crucial parameter for the adequate description of the water balance by the models.

INTRODUCTION

During the last decades a significant increase of water bodies contaminated with pesticides has been observed for various areas worldwide. Due to the use of aquifers for drinking water supply, an assessment and a sound prediction of the impact of agricultural practice concerning pesticide application is mandatory. In the framework of the EU-Project PEGASE (Pesticides in European Groundwaters: detailed study of representative aquifers and simulation of possible evolution scenarios) the test area 'Zwischenscholle' (20 km²), located 30 km west of Cologne (Germany), was selected for a joint modelling on pesticide transport at regional scale with four different modeling approaches. Two of these four models are MARTHE and TRACE, which calculate water flow in variable saturated porous media using a generalized Richards' equation. The other two models, MACRO and ANSWERS are used in a coupled approach, linking both aforementioned unsaturated zone models with a groundwater flow model. The validation of water flow and transport processes at regional scale is rather difficult. Thus in a first step a lysimeter data set was chosen to check the modelling of water flow and pesticide transport at local scale. Cropped lysimeters were used to monitor water flow and fate of the pesticide methabenzthiazuron (MBT) for a period of almost two years.

The lysimeter station was located in the test area 'Zwischenscholle' and the lysimeter soil approximates some of the soils found in the test area. The experimental data of the lysimeter were supposed to validate the modelling of plant related processes, water flow and pesticide transport in the unsaturated zone at local scale.

Large undisturbed lysimeters are a common experimental setup for investigations of pesticide transport [START_REF] Bergström | Use of lysmeters to estimate leaching of pesticides in agricultural soils[END_REF][START_REF] Boesten | Simulation of bentazon leaching in sandy soil from Mellby (Sweden) with the PESTLA model[END_REF][START_REF] Keller | Mobility and Dissipitation of 14C-labeled atrazine, Metolachlor, and Primisulfuron in undisturbed field lysimeters of a coastal plain soil[END_REF][START_REF] Vink | Simulation and model comparison of unsaturated movement of pesticides from a large clay lysimeter[END_REF][START_REF] Schoen | Preferential flow and solute transport in a large lysimeter, under controlled boundary conditions[END_REF][START_REF] Mikata | Mobility and degradation of the herbicide imazosulfuron in lysimeters under field conditions[END_REF], they are particularly applied for pesticide registration purposes.

The main advantages of lysimeters are the controlled boundary conditions and the measurement of actual evapotranspiration, soil moisture and drainage for a large soil volume.

In addition the pesticide concentrations of the drainage water can be determined. The representativeness of lysimeter observations for the field scale behaviour of transport processes is still debated. Compared with field measurements with suction plates [START_REF] Jene | Comparison of lysimeters to the actual field situation with regard to the movement of soil water and a 14 C-labelled herbicide and bromide. The environmental fate of Xenobiotics[END_REF] measured 40% more bromide leaching in lysimeters. Comparative modeling of pesticide transport has been the subject of several studies. Model comparisons on pesticide transport have been carried out for lysimeters (Bergtröm and Jarvis, 1994;[START_REF] Vink | Simulation and model comparison of unsaturated movement of pesticides from a large clay lysimeter[END_REF][START_REF] Francaviglia | Comparing and evaluating pesticide leaching models: results for the Tor Mancina data set (Italy)[END_REF] as well as for field studies [START_REF] Pennel | Evaluation of five simulation models for predicting aldicarb and bromide under field conditions[END_REF]Diekkrüger et al. 1995a;[START_REF] Armstrong | Comparison of pesticide leaching models: results using the Brimstone Farm data set[END_REF][START_REF] Gottesbüren | Comparison of pesticide leaching models: results using the Weiherbach data set[END_REF][START_REF] Tiktak | Application of pesticide simulation models to the Vredepeel dataset II Pesticide fate[END_REF]Vanclooster and Boesten, 2000). The most recent and probably most extensive model comparsion was summarized by Vanclooster et al. (2000). The major outcome of this comparative pesticide modelling on lysimeter and field data by [START_REF] Armstrong | Comparison of pesticide leaching models: results using the Brimstone Farm data set[END_REF], [START_REF] Francaviglia | Comparing and evaluating pesticide leaching models: results for the Tor Mancina data set (Italy)[END_REF], [START_REF] Gottesbüren | Comparison of pesticide leaching models: results using the Weiherbach data set[END_REF], [START_REF] Tiktak | Application of pesticide simulation models to the Vredepeel dataset II Pesticide fate[END_REF] and Vancloster and Boesten (2000) is what Diekkrüger et al. (1995a) also stressed: The influence of the modellers' experience on model results is large, probably larger than the influence of the selected model concept. [START_REF] Bergström | Evaluation and comparison of pesticide leaching models for registration purposes[END_REF] found very similar results for the five models included in their comparison.

They also note that, besides from taking the relevant processes into account, the identification of correct model input parameters plays a key role in predicting pesticide transport.

Methabenzthiazuron is an effective herbicide in grain and certain vegetable crops. It is classified as a rather persistent compound [START_REF] Rouchaud | Metabolism of methabenzthiazuron in the soil of pea crops[END_REF]) with a high sorptivity (Diekkrüger et al., 1995b;[START_REF] Wüstemeyer | Abbau und Transport von Ethidimuron und Methabenzthiazuron in einer Parabraunerde-Beitrag zur Validierung des Lysimeterkonzepts[END_REF].

MARTHE and TRACE are rather recent developments, based on the extension of codes originally developed for the description of groundwater transport at regional scales. TRACE comes from 3DFEMWATER (Yeh at al., 1992), which allowed to calculate both saturated and unsaturated water flow. MACRO is a classical 1-dimensional model. Since the introduction of MACRO it has already been used in several model comparison studies on pesticide fate modelling [START_REF] Bergström | Evaluation and comparison of pesticide leaching models for registration purposes[END_REF]Diekkrüger et al. 1995a;[START_REF] Vink | Simulation and model comparison of unsaturated movement of pesticides from a large clay lysimeter[END_REF]Vanclooster et al., 2000), where it has proven a broad applicability to several pesticide transport problems. During this study the version 5.0 is applied. ANSWERS is a capacity based regional scale model for the vadose zone transport of solutes. It has already been evaluated in the context of agricultural non-point source water quality models [START_REF] Kosky | Evaluation of three distributed parameter hydrologic/water quality models[END_REF]. The model concepts are rather contrasting (Table 1). MACRO is the only model accounting for preferential flow and transport. Except for ANSWERS all models are based on the Richards' equation for the calculation of soil water flow. MACRO and ANSWERS solve the convection-dispersion equation (CDE) with a common numeric scheme, whereas TRACE is coupled with 3dLEWASTE, which applies a hybrid Lagrangian/Eulerian method to solve the CDE and MARTHE uses the total variation diminishing method. For pesticide degradation only TRACE coupled with 3dLEWASTE uses a simple first-order kinetics approach, whereas all the other models allow the use of a soil temperature and moisture dependent biodegradation approach.

Basically, there are three main objectives: (i) to identify the crucial plant parameters and test the applicability of the relevant plant data base, (ii) to check the model performances on water flow and MBT transport and (iii) to compare the performance of the different model concepts.

METHODS AND DATA

Experimental data

Five undisturbed soil monoliths (free draining lysimeters) containing an Orthic Luvisol were used to monitor the soil water balance. The sampling distance between the monoliths was roughly 1 m. The monoliths were 1.1 m long and had a surface area of 1.0 m 2 . Three pedogenetic soil horizons were distinguished (Table 2). One of the lysimeters was treated with a dose of 2.8 kg ha -1 as TRIBUNIL  corresponding with an active ingredient application of 248.11 mg m -2 [phenyl-U-14 C]methabenzthiazuron during the pre-emergence period of winter wheat at the 25 th of November 1988. The winter wheat (Triticum Aestivum L.) was harvested and soil samples were taken with a hand auger at 3 rd of August 1989 (252 days after  Reg. Trademark, Bayer AG application). In the following vegetation period winter barley (Hordeum Vulgare L.) was cropped and harvested on the 11 th of May 1990. The next vegetation was oat, which was harvested at the 13 th of August 1990, when a last soil sampling was carried out 627 days after the application.

From the 25 th of November 1988 until the 13 th of August 1990 the meteorological parameters precipitation, air humidity, air temperature, wind speed and radiation were monitored on a daily basis (Figure 1). Roughly every three days soil moisture was measured with a neutron probe. Leachate was collected on a three weekly basis. The amount of lysimeter drainage can be measured directly while the actual evapotranspiration Et a (mm) was calculated from the soil water balance,

     D P Et a [1]
where P is the precipitation (mm), D is the drainage (mm) and Δθ is the change of soil moisture (mm) during the considered period of time.

MBT was extracted from the soil with acetone/ethyl acetate/chloroform. Using thin-layer chromatography the detection limit of 14 C-labeled MBT was 0.67 µg l -1 .

TRACE/3dLEWASTE

TRACE [START_REF] Vereecken | TRACE A mathematical model for reactive transport in 3D variably saturated porous media[END_REF] calculates the 3-dimensional unsaturated/saturated water flow in porous media. A modified Picard-iteration scheme [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF] is applied in combination with a preconditioned conjugate gradient method in order to solve the following generalized 3-dimensional Richards` equation numerically,

) , ( )) , ( ) , ( ( ) , ( t s Q h s H h s K t t s h F h           with h h S h F      [2]
where K is the hydraulic conductivity (LT -1 ),  s is the position vector in a three dimensional space, h is the pressure head (L), H is the total head (L), Q is the source/sink term (T -1 ), θ is the soil water content (L 3 L -3 ), t is time (T), S h is the specific storage coefficient (L -1 ) and F h is the specific water capacity (L -1 ). For the spatial discretization hexahedral Galerkin-type finite elements are used. In order to take plant related processes into account the crop growth module SUCROS (Simplified and Universal Crop growth Simulation, [START_REF] Spitters | A simple but universal crop growth simulation model[END_REF] is implemented. In contrast to many other plant modules SUCROS estimates the assimilation rate from plant specific photosynthesis parameters and radiation. The calculated increase in biomass is used to predict leaf area growth and leaf area index (LAI), which is therefore not model input. The required plant data values can e.g. be found in [START_REF] Van Heemst | Plant data values required for simple crop growth simulation models: review and bibliography[END_REF]. The crop coefficients (K c -values) for the scaling of the reference evapotranspiration can be assigned for different seasonal stages according to the approach of Doorenbos and Pruitt (1978). Based on LAI the potential evapotranspiration is split into the fractions of potential evaporation and potential transpiration according to Beer's law. The actual transpiration is calculated from the potential transpiration in dependence of the pressure head according to the approach of [START_REF] Feddes | Simulation of field water use and crop yield[END_REF]. For the soil evaporation a flux boundary is applied to the uppermost element until the given pressure head h min is reached (e.g. -10 -4 cm). At this point TRACE switches to a fixed head boundary condition set to h min .

3dLEWASTE (Yeh et al., 1992) is a hybrid Lagrangian-Eulerian finite element model of reactive solute transport through unsaturated/saturated media. 3dLEWASTE numerically solves the Lagrangian form of the convection-dispersion transport equation (CDE),

C D C V Dt DC          ) ( [3]
where C is the solute concentration (M L -3 ), DC/Dt is the material derivative of C with respect to time t,  D is the dispersion coefficient tensor (L 2 T -1 ) and  V is the pore velocity vector (L T -1 ) for the x,y and z direction. Here the advective term is solved in a mobile (Lagrangian) coordinate system using a single step reverse particle tracking, while the diffuse term is solved in a fixed (Eulerian) coordinate system. A backward differencing scheme in time is applied. A more detailed description of the hybrid Lagrangian-Eulerian approach can be found at Yeh et al. 1992. For microbial decay a first order degradation rate coefficient is applied to the sorbed and the liquid concentrations. A linear, Freundlich or Langmuir isotherm can be applied for sorption. For the coupling between TRACE and 3dLEWASTE a file is used containing the Darcy fluxes and the water contents for every finite element at every time step. For simplicity the coupled model TRACE and 3dLEWASTE is referred to as 'TRACE' during the following text.

ANSWERS

The Areal Non-point Source Watershed Environmental Response Simulation, ANSWERS [START_REF] Bouraoui | A distributed physical approach for surface-subsurface water transport modeling in agricultural watersheds[END_REF], is originally a watershed scale, diffuse pollution model for long term simulation. The core of the system is a one-dimensional vertical model based on a capacity approach for the soil water flux. A variable vertical segmentation is considered to account for water movement through the soil profile. Infiltration into the uppermost layer is simulated by the Green and Ampt equation [START_REF] Green | Studies in soil physics. I. The flow of air and water through soils[END_REF]. Soil water redistribution from upper layer of soil to the root zone (discretized in 9 layers) is determined with the use of a Brooks and Corey type equation [START_REF] Brooks | Hydraulic Properties of Porous Media[END_REF] on the basis of vertical downwards gravity flow, with a hydraulic conductivity related to the average water content of the upper layer.

Similarly, percolation from root zone to the underlying unsaturated zone (drainage) is determined with a hydraulic conductivity related to the average water content of the root zone.

The main parameters describing the hydraulic properties of the soil are the saturated hydraulic conductivity K s (cm h -1 ), the porosity Φ (cm 3 cm -3 ), the residual water content  r (cm 3 cm -3 ), the pore size distribution index λ b (-), the bubbling pressure ψ b (cm) and the field capacity F c (cm). They are obtained from soil texture and organic matter content with the Pedotransfer functions of [START_REF] Rawls | Estimation of soil retention and hydraulic properties[END_REF].

Soil evaporation and plant transpiration are modelled separately using Ritchie's equation [START_REF] Ritchie | Model for predicting evaporation from a row crop with incomplete cover[END_REF], where soil evaporation is related to the soil moisture of the upper layer and the leaf area index LAI (m 2 m -2 ). Plant transpiration is extracted from root zone assuming a uniform root profile. The parameters describing the plant behaviour in terms of water uptake and actual evapotranspiration are obtained from a database including 78 different types of crop [START_REF] Knisel | GLEAMS Groundwater loading effects of agricultural management systems[END_REF]. From sowing to harvest values of the LAI and root depth are given for phenological stages.

The transport of solutes is calculated with the convection-dispersion equation. For sorption, degradation and plant root uptake of pesticides the approach of ANSWERS is similar to GLEAMS [START_REF] Knisel | GLEAMS Groundwater loading effects of agricultural management systems[END_REF]. Two important assumptions are made: the degradation is only a function of temperature (zero below 0°C and maximum at 25°C) and plant uptake is 10% of water uptake. For the degradation of pesticides alternatively the equation of [START_REF] Graham-Bryce | Movement of fluometuron and 36 Clin soil under field conditions: Computer simulations[END_REF] accounting for soil temperature and soil moisture dependent decay (Eq. 9) can be used.

MACRO

MACRO 5.0 [START_REF] Jarvis | Macro5.0 XXXXXXX[END_REF]) is a one-dimensional dual-permeability model, operating at the scale of a soil profile. The model accounts for preferential flow and transport in soil macropores by dividing the soil pore system into two parts, one part with a high flow and low storage capacity (macropores) and the remainder with a low flow and a high storage capacity (micropores). The boundary between the pore regions is defined by the fixed pressure head C ten (L) having a corresponding water content θ b (L 3 L -3 ) and corresponding hydraulic conductivity K mi (L T -1 ). The 1-dimensional form of the Richards equation is used to model flow in the micropores,

) , ( 1 ) , ( t z Q z h h z K z h h F h                        with h F h     [4]
where z is the vertical coordinate (L). Eq. 4 is solved for finite differences using the implicit iterative procedure proposed by [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF]. Soil water retention and unsaturated hydraulic conductivity are calculated using a modified form of the Mualem/van Genuchten approach [START_REF] Van Genuchten | A closed form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF][START_REF] Vogel | Effect of the shape of the soil hydrulic functions near saturation on variably-saturazed flow predictions[END_REF] accounting for the macropore/micropore dichotomy by using the boundary pressure head partitioning the total porosity into micropores and macropores [START_REF] Wilson | Modeling the hydraulic properties of a multiregion soil[END_REF][START_REF] Mohanty | New piecewise-continuous hydraulic functions for modelling preferential flow in an intermittent flood-irrigated field[END_REF]. Flow in the macropores is calculated using the kinematic wave equation [START_REF] Germann | Kinematic wave approach to infiltration and drainage into and from soil macropores[END_REF], assuming gravity-dominated flow (i.e. neglecting capillarity). The hydraulic conductivity function in the macropores is given as a simple power law expression of the macropore degree of saturation.

The one-dimensional convection-dispersion equation is applied for solute transport in the micropores,

                       i m b ne mac m U qC z C D z t S f f C    1 [5]
where the source/sink term U i (M L -3 T -1 ) represents different processes like mass exchange between flow domains, kinetic sorption, solute uptake by the crop, biodegradation and lateral leaching losses to drains or groundwater. S is the sorbed concentration (M L -3 ), f mac is the mass fraction of solid material in contact with water in the macropores (-), f ne is the fraction of the solid material providing non-equilibrium sorption (-), ρ b is the soil bulk density (M L -3 ), θ m is the mobile water content (L 3 L -3 ), q is the water flow (L T -1 ) and D is the Dispersion coefficient (L 2 T -1 ). Transport in macropores is calculated neglecting dispersion-diffusion, but accounting for adsorption by the parameter f mac (-) that partitions the sorption constant between the two flow regions. Diffusive mass exchange between the two pore regions is calculated using approximate first-order equations based on an effective diffusion path length a scale (L). Solute transport is solved by a Crank-Nicolson finite difference scheme utilizing an iterative, fully upstream weighting procedure with an empirical correction for numerical dispersion.

Root water uptake is calculated from evaporative demand, root distribution and soil moisture using a modified version of the approach developed by [START_REF] Feddes | Simulation of field water use and crop yield[END_REF] accounting for water stress compensation [START_REF] Jarvis | A simple empirical model of root water uptake[END_REF]. It is assumed, the crop can adjust to stress in one part of the root system by increasing uptake from other parts where the soil moisture conditions are more favourable. Root density is assumed to be distributed logarithmically with depth. Beer`s law is used to partition the potential evapotranspiration into one fraction transpirated by the canopy and the remaining fraction of evaporation from the soil. This is based on the green and on the total leaf area indices, given as a function of the day number in the year as user-specified input.

The heat conduction equation is solved using a standard Crank-Nicholson finite difference scheme. The effect of soil moisture and temperature on the first-order kinetics degradation of pesticides is estimated with the approach of Boesten and van der Linden (1991):

  ref T T ref ref e                 for T > 5°C [6]
where λ is the degradation rate coefficient (d -1 ), λ ref is the reference rate coefficient (d -1 ), θ ref is the reference water content (cm cm -3 ), T is the soil temperature (°C), T ref is the reference temperature (°C), and α (°C -1 ) and β (-) are empirical parameters. The pesticide uptake by roots is modelled as a function of root water uptake and pesticide concentration. An empirical concentration factor is used to define the fraction of pesticide concentration taken up by the roots.

MARTHE

MARTHE (Modelling Aquifers with Rectangular cells, Transport and hydrodynamics) was

originally developed as a 3-dimensional groundwater model designed to compute water flow and solute transport in saturated porous media [START_REF] Thiéry | Modélisation 3D du transport de masse avec le logiciel MARTHE, version 5.4[END_REF]. Additional routines allow the computation of unsaturated flow, thus MARTHE also solves the 3-dimensional form of the Richards equation numerically. The spatial discretization is based on finite differences. The resulting matrix equation is iteratively solved using conjugate gradients combined with Choleski pre-conditioning. Advective, diffusive and dispersive transport can be simulated using three different techniques: CDE based on finite differences, Total Variation Diminishing and Method of Characteristics.

The main features of the plant module are very similar to MACRO. The evolution of the LAI is a function of main development stages (i.e. germination, maturity and harvest). Like for MACRO an exponential function is used to divide the potential evapotranspiration into transpiration and evaporation. Canopy interception is not considered. For the root system development of crops a linear growth from germination until maturity is assumed. Between maturity and harvest the specified maximum root depth remains constant. The root density can be calculated from several functions. During this study an exponentially decreasing root density was assumed. For the estimation of the actual transpiration the reduction according to the water stress compensation concept [START_REF] Jarvis | A simple empirical model of root water uptake[END_REF]) is taken into account.

The root uptake of solutes is simply calculated from a solution uptake factor (-) describing the fraction of mass lost to root uptake by plants. Concerning solute degradation a sequential first-order decay, a temperature and soil moisture dependent first-order decay or a simple first-order decay can be taken into account. For the soil temperature and moisture dependent first-order decay a parameterized Graham-Bryce approach [START_REF] Graham-Bryce | Movement of fluometuron and 36 Clin soil under field conditions: Computer simulations[END_REF] or the concept of Boesten and van der Linden (1991) can be used.

Model input

The functional relation between pressure head, soil water content and unsaturated hydraulic conductivity plays a key role for the modeling of water flow [START_REF] Vereecken | Analysis of multi-step outflow data and pedotransfer functions to characterise soil water and solute transport at different scales[END_REF][START_REF] Herbst | The influence of the spatial structure of soil properties on water balance modeling in a microscale catchment[END_REF]. For MARTHE and TRACE the soil water retention function of van Genuchten (1980) with parameter m=1 was used, which is equivalent to the equation proposed by [START_REF] Brutsaert | Probability laws for poresize distributions[END_REF],

n r s r h h ) ( 1 ) (          [7]
where  denotes the water content (cm 3 cm -3 ), h is the pressure head (cm),  r is the residual water content (cm 3 cm -3 ),  s is the water content at saturation (cm 3 cm -3 ),  is the inverse of the bubbling pressure (cm -1 ) and the dimensionless shape parameter n (-).

By using m=1 instead of m=1-1/n the closed analytical expression of the Mualem/van Genuchten approach (van Genuchten, 1980) for the K(h) function is lost. Therefore the unsaturated hydraulic conductivity function of [START_REF] Gardner | Some steady state solutions of the unsaturated moisture flow equation with application to evaporation from a water table[END_REF] was applied,

c s h b K h K ) ( 1 ) (   [8]
where K is the unsaturated hydraulic conductivity (cm d -1 ), K s is the saturated hydraulic conductivity (cm d -1 ), and b (cm -1 ) and c (-) are empirical parameters. The soil hydraulic properties listed in Table 4 were determined by [START_REF] Vereecken | Analysis of multi-step outflow data and pedotransfer functions to characterise soil water and solute transport at different scales[END_REF]. They derived the soil hydraulic properties from the soil properties listed in Table 2 with the Pedotransfer functions (PTF) of [START_REF] Vereecken | Estimating the soil moisture retention characteristic from texture, bulk density and carbon content[END_REF][START_REF] Bergström | Use of lysmeters to estimate leaching of pesticides in agricultural soils[END_REF] and fitted the  r and  values. TRACE and MARTHE allow for the use of the modified Mualem/van Genuchten approach described above, while MACRO is based on the common Mualem/van Genuchten approach (with m=1-1/n) modified for the dual porosity approach [START_REF] Wilson | Modeling the hydraulic properties of a multiregion soil[END_REF][START_REF] Mohanty | New piecewise-continuous hydraulic functions for modelling preferential flow in an intermittent flood-irrigated field[END_REF][START_REF] Vogel | Effect of the shape of the soil hydrulic functions near saturation on variably-saturazed flow predictions[END_REF]. In order to obtain comparable soil hydraulic functions for the Richards' equation based models the retention and unsaturated hydraulic conductivity functions for MACRO (Mualem/van Genuchten, m=1-1/n) were derived in a two step procedure. First  r ,  s , , n and K s were fitted to the already available functions of [START_REF] Vereecken | Analysis of multi-step outflow data and pedotransfer functions to characterise soil water and solute transport at different scales[END_REF] using the least squares procedure of RETC [START_REF] Van Genuchten | The RETC Code for Quantifying the Hydraulic Functions of Unsaturated Soils[END_REF]. In a second step the parameters defining the macropore system the saturated matrix conductivity K mi (cm d -1 ), the diffusion path length a scale (cm) and the boundary pressure head C ten (cm) were calibrated on the measured MBT concentrations in the drainage water (Table 4). For ANSWERS soil hydraulic parameters like K s and field capacity were calculated internally with the PTF of [START_REF] Rawls | Estimation of soil retention and hydraulic properties[END_REF] from the properties given in Table 2. The use of other soil hydraulic functions is not very appropriate, since the capacity based approach implemented in ANSWERS is closely linked to the PTFs of [START_REF] Rawls | Estimation of soil retention and hydraulic properties[END_REF]. The initial values of soil moisture or pressure head for the model start were derived from the neutron probe measurements.

All the models can handle a linear, Freundlich or Langmuir sorption isotherm. For all of the models the sorption of MBT was described with a linear sorption isotherm. The values of the distribution coefficient K d (cm 3 g -1 ) are calculated from a partition coefficient K oc of 527 cm 3 g -1 , found in the Agritox database (http://www.inra.fr/agritox), and the organic matter content assuming K d =K oc f oc . The resulting K d for the three soil horizons from top to the bottom is 5.27, 2.11 and 1.58 cm 3 g -1 respectively. Thus the retardation factor for the top horizon at water saturation is 21.2, indicating a high sorption of MBT. For MARTHE and ANSWERS a temperature and soil moisture dependent degradation based on the equation of Graham-Bryce

et al. (1982) is taken into account, )] ( ) ln( [ 1 ) , (       T g M b a e M T  [9]
where λ is the degradation rate (d -1 ), T is the soil temperature (°K), M is the gravimetric soil moisture (%) and a (-), b (-) and g (°K) are parameters. From batch experiments [START_REF] Wüstemeyer | Abbau und Transport von Ethidimuron und Methabenzthiazuron in einer Parabraunerde-Beitrag zur Validierung des Lysimeterkonzepts[END_REF] found the following parameters related to MBT:

1 4476 ) ln( 343 . 1 0 . 7 ) ln(       T M  [10]
From which the half-life (DT 50 ) can be calculated as follows: and θ ref =0.39) of 162 days was calculated according to the approach of [START_REF] Boesten | Modelling the influence of sorption and transformation on pesticide leaching and persistence[END_REF]. This DT 50 is applied for MACRO. Figure 2 shows a comparison of the Graham-Bryce approach with the Boesten and van der Linden approach. It is clearly visible, that the Boesten and van der Linden approach with a reference DT 50 of 162 days gives higher half-lifes than the Graham-Bryce approach. Thus the remaining mass of MBT in soil should be smaller for MARTHE and ANSWERS applying the Graham-Bryce approach than for MACRO and TRACE, where the DT 50 is generally higher. For all the models the degradation rate for sorbed and dissolved MBT is supposed to be the same. The dispersion length was assumed to be 1.7 cm for all of the models.

MACRO, MARTHE and ANSWERS calculate the soil temperature from daily mean air temperature by solving the equation of heat diffusion and convection. For MARTHE the mineral thermal conductivity is set to 1.5 W m -1 °C-1 , water thermal conductivity is set to 0.6 W m -1 °C-1 , mineral volumetric specific heat is set to 2•10 6 J m -3 °C-1 and water specific heat is set to 4185 J kg -1 °C-1 . MACRO applies the approach of [START_REF] Jansson | Simulation model for soil water and heat conditions. Description of the SOIL model[END_REF] to estimate thermal conductivity and volumetric heat capacity from the soil properties summarized in Table 2.

The potential reference evapotranspiration ET p (Figure 1) was calculated according to the approach of Penman/Monteith [START_REF] Monteith | Vegetation and the atmosphere[END_REF]. Potential evapotranspiration and precipitation are applied as the upper boundary condition, whereas a seepage face is applied to the lower boundary at the bottom of the lysimeter. The seepage face boundary is characterized by a no-flow boundary for unsaturated conditions:

0 ) , (  t z q for 0  h [12]
If the seepage face becomes saturated, the boundary turns into a prescribed head boundary with h(z,t)=0. The vertical lysimeter wall is set as a no-flow boundary condition.

For MACRO, MARTHE and TRACE the spatial discretization consists of 110 elements of 1 cm thickness each. The 3-d models MARTHE and TRACE are run in a 1-d mode. Due to the infiltration approach according to [START_REF] Green | Studies in soil physics. I. The flow of air and water through soils[END_REF] the spatial discretization for ANSWERS requires a thickness of 40 cm for the uppermost element. The eight elements beneath have a thickness of 8.75 cm each.

Validation criteria

A commonly used criterion for model validation is the root mean square error (RMSE), where the root of the mean squared residuals is calculated. The RMSE has the unit of the considered variable. The squared residuals are also used for the second criterion applied, which is the coefficient of model efficiency CME [START_REF] Nash | River flow forecasting through conceptual models: Part I -A discussion of principles[END_REF]. Here they are used to determine the proportion of the deviation from the observed mean, which can be explained by the model,

           n i i omean o n i i s o n i i omean o x t x t x t x x t x CME 1 2 1 2 1 2 ) ) ( ( )) ( ) ( ( ) ) ( ( [13]
where x o is the observed value, x s is the simulation result at time t and x omean is the arithmetic mean of the observed values. The CME is a dimensionless criterion that can have negative values. The highest value possible is 1, indicating that observation and model are completely in agreement. The Index of Agreement IA [START_REF] Willmott | On the validation of models[END_REF] is also dimensionless and ranges between 0 and 1:

                      n i i omean s omean o n i i s o x t x x t x t x t x IA 1 2 1 2 1 [14]
Because CME and IA are dimensionless, they can be used to compare the model quality between different variables, while the RMSE gives an idea about the model error in the units of the variable under consideration.

RESULTS

Water flow

One aim of this study was to identify the crucial parameters of plant related processes concerning the pesticide transport in soil. During the calibration process of the four models three crucial parameters were identified: The leaf area index LAI, the root depth Rd and the crop conversion factor K c . Table 3 summarizes the calibrated and the original maximum values of this three parameters. Only for MACRO the maximum LAI was slightly modified during calibration. For all the crops the emergence and harvest dates were known, thus the uncertainty concerning the LAI development was rather small and modifications of the LAI did not improve model results significantly. A quite large variability of the maximum LAI was detected for the different plant databases. The difference between MACRO and MARTHE, which are based on the same plant database, and ANWSERS is larger than the difference between the three crops. Compared to the plant database LAIs the LAI estimated with TRACE for oat is much smaller. The LAI estimated for winter wheat and winter barley are much more similar. The root depth defines the zone of the soil, where the sink term of transpiration influences the soil moisture. Differences between calibrated and uncalibrated maximum Rd range between 5 and 40 cm, for MACRO and ANSWERS respectively. The crop conversion factor K c determines the amount of ET p . Among this three plant parameters the K c has probably the largest influence on the model results and the largest differences between calibrated and uncalibrated plant parameters can be detected for the this parameter.

The variation of the calibrated maximum K c for different crops is large (Table 3). The lowest K c was found for the winter barley, while the highest values were assigned to winter wheat.

With ANSWERS the LAI and K c factors were unchanged during calibration. For TRACE and MACRO the Rd and K c -values were significantly modified during the calibration procedure.

The same basically holds true for MARTHE, but the K c was changed more moderately. For MARTHE a constant K c was applied for every crop, without any temporal variability according to phonological stages.

The total water balance (Table 6) is well reproduced by all models. ANSWERS slightly overestimates the total amount of actual evapotranspiration ET a , whereas the other models slightly underestimate the total ET a . The relative errors for single water balance components are small, except for the change in soil moisture, which is merely a result of the small absolute amount of change in soil moisture.

Apart from the overall water balance, the reproduction of the temporal evolution of the water balance components is relevant. Figure 3 shows a comparison between observed and modelled cumulative actual evapotranspiration. In general the four models match the measurements, which exhibit only small standard deviations, except for the two drying periods with high evapotranspiration demands. During the first spring period (1989) all the models applied slightly underestimate the amount of actual evapotranspiration. This is vice versa for the second period (spring 1990), when MARTHE shows an increase in actual evapotranspiration too early and too low. ANSWERS reproduces this increase too early and in total too high. The calculation of the coefficient of model efficiency (CME) with the mean of the measured actual evapotranspiration and the corresponding model results reveals that ANSWERS, MACRO and MARTHE are close to each other in their ability to reproduce evapotranspiration while TRACE shows the highest CME (Table 5). For the models the ratio between cumulative actual evapotranspiration (ET a ) and cumulative potential evapotranspiration (ET p ) differs for several periods. The K c -values applied already indicate this. Combined with the different reduction methods to estimate ET a from ET p for the different models the ratio between ET a and ET p shows large deviations. For example at the end of the first drying period during summer 1989 at the same day the ratio ET a /ET p ranges between 0.65 and 1.0, for TRACE and MARTHE respectively, whereas the resulting ET a is quite similar for all models (Figure 3).

A comparison of calculated and observed volumetric soil moisture at two depths is shown in

Transport of Methabenzthiazuron

In the models MARTHE, MACRO and ANSWERS degradation is soil moisture and temperature dependent. Figure 6 shows a comparison between the model results concerning soil temperature at 10 cm depth. This depth was chosen because the measured soil temperature at 10 cm depth exhibits a high amplitude, which is a result of being close to the soil surface, where the soil heat flux is mainly driven by the atmospheric conditions.

Furthermore most of the degradation takes place in the uppermost soil layer with high organic matter content. The modelled soil temperatures of MARTHE, MACRO and ANSWERS are very close to each other and they are close to the measurements. The temporal variability of soil temperature is well reproduced by the models, which is supported by the small RMSE (Table 5). For MARTHE, MACRO and ANSWERS the degradation was calculated with the modelled soil temperature and soil moisture. Figure 7a reveals that 252 days after application the remaining total MBT mass was found in the upper 10 cm. This observation is generally reproduced with TRACE, MACRO and MARTHE, although the three models show a small amount of MBT in the soil layer between 10 and 20 cm depth. The amount of MBT in the upper 10 cm estimated with MARTHE is very close to the measurement while TRACE clearly overestimates the remaining mass of MBT. Related to the applied mass, TRACE estimated 12 % too much. MACRO also clearly overestimates the MBT residues for the first sampling date. With ANSWERS no MBT is found in the upper compartment (0-40 cm depth).

The total mass of MBT is found in the depth between 40 and 70 cm. Related to the applied mass, the total mass left in the profile according to ANSWERS is 13.6 % too much.

After 627 days the measurements show that the total mass left in the profile is just little lower than after 252 days and it is divided into the two uppermost layers (0-20 cm). Compared to the measurements the results obtained with MACRO and TRACE are the closest (Figure 7b).

MACRO and TRACE underestimate the amount of MBT, but they reproduce the right depth.

MARTHE also reproduces the right depth, but the degradation of MBT is clearly overestimated. The same holds for ANSWERS and additionally the mass is again estimated deeper in the profile than measured.

As assumed from Figure 2 the estimated mass of MBT residuals is higher for both sampling dates for MACRO and TRACE than for MARTHE and ANSWERS. This is an effect of the higher DT 50 . For TRACE the DT 50 is constantly 200 days, while MACRO modifies the DT 50 of 162 days according to soil moisture and temperature following the approach of Boesten and van der Linden (1991). MARTHE and ANSWERS compute smaller DT 50 from the Graham-Bryce equations (Eq. 10 and 11). For the first sampling date 252 days after application the RMSE calculated for the 11 depths (see Figure 7) is 30.8, 3.9, 9.7 and 11.7, for ANWERS, MARTHE, TRACE and MACRO respectively. For the second sampling 627 days after application the RMSE for the 11 depths is 13.7, 10.1, 5.1 and 6.4, respectively. The models applied with a higher DT 50 (TRACE and MACRO) reproduce the measurements of the second date better than the ones of the first date. This is vice versa for MARTHE, which is based on the Graham-Bryce approach. MARTHE reproduced the measurements of the first sampling better than for the second sampling.

According to Figure 8 the measured concentrations of MBT in the drainage water show a small peak roughly 100 days after application. The total mass of MBT lost by leaching during the experimental period is 14.6 µg m -2 , equivalent to 0.0059 % of the applied mass of MBT.

None of the models predicts MBT in drainage water, except the calibrated MACRO estimates at total amount of MBT leaching of 7.8 µg m -2 . MACRO was mainly calibrated on the macropore parameters listed in Table 4. Figure 8 shows that the peak of MBT leaching predicted with MACRO is a little delayed, and the total amount of MBT leaching is underestimated. Apart from the calibration of the macropore parameters with the water balance components the fraction of sorption sites in the macropores was calibrated with the measured MBT leaching, f mac was set to 0.005 (-), which corresponds to the small macropore volume fraction (θ s -θ b , Table 4).

DISCUSSION

Three crucial parameters were identified during the calibration procedure. The calibration of the LAI was moderate, probably due to the exactly known emergence and harvest dates. A striking point is the higher variability between different plant databases than between calibrated and uncalibrated model parameters. A query at the plant parameter database (PlaPaDa, [START_REF] Breuer | PlaPaDaan online plant parameter data drill for eco-hydrological modelling approaches[END_REF] revealed LAIs between 2.3 and 4.6 for winter wheat. This

shows the variability of this parameter for different environmental conditions, and it might indicate that the model functions are adapted to the related plant database. Like the LAI the root depths were rather moderately calibrated, except for the root depth of oat. For this crop three of the four modellers increased the root depth. From the calibration of LAI and root depth it could be assumed that the chosen databases of [START_REF] Knisel | GLEAMS Groundwater loading effects of agricultural management systems[END_REF], [START_REF] Van Heemst | Plant data values required for simple crop growth simulation models: review and bibliography[END_REF] and [START_REF] Focus | FOCUS groundwater scenarios in the EU review of active substances[END_REF] in combination with the chosen model are generally applicable to the Zwischenscholle test site. The really important plant parameter is the crop conversion factor K c . This parameter has a strong influence on the overall water balance. Except for ANSWERS all models were calibrated on K c for each of the crops. For winter wheat and oat the K c values were increased, while for winter barley K c was decreased. The decrease of the K c for winter barley was required to account for the effect of a plant disease keeping the barley from a normal plant development. Against this background and having in mind the amount of calibration of K c even for this parameter the selected databases are basically applicable for the Zwischenscholle. But care should be taken not to underestimate the actual evapotranspiration by using K c values that are too small. This can only be checked in terms of the water balance.

Further there might be, compared to the field situation, generally a slightly higher evapotranspiration for lysimeters. [START_REF] Bergström | Evaluation and comparison of pesticide leaching models for registration purposes[END_REF] and [START_REF] Boesten | Simulation of bentazon leaching in sandy soil from Mellby (Sweden) with the PESTLA model[END_REF] attribute this to the 'oasis effect', i.e. the lysimeters are partially surrounded by hard surfaces without any evaporation, which presumably causes in combination with wind a lateral flow of dry air over the lysimeters, thereby increasing the evapotranspiration of the lysimeters.

As expected the soil moisture near the soil surface is clearly influenced by evapotranspiration. This is reflected in the model results. If the model correctly reproduces the amount of transpiration and evaporation the soil moisture at the depth of 25cm is well described.

Supported by the CME and IA all the model results for evapotranspiration and soil moisture at 25 cm depth are clearly acceptable, which is basically the result of well described plant related processes from the calibration of plant parameters. The soil moisture at the depth of 85cm is influenced by root water uptake and drainage. Even for the soil moisture at 85 cm a proper reproduction can be stated for all the models (Table 5), whereas the temporal course of drainage is reproduced with different quality. According to [START_REF] Nash | River flow forecasting through conceptual models: Part I -A discussion of principles[END_REF] According to the validation criteria TRACE reproduces the soil moisture at 85cm even slightly better than MARTHE. But the drainage is described significantly better by MARTHE. The validation criteria for the soil moisture at 85cm calculated with ANSWERS are in the same range for all the models, while the drainage is not that well reproduced. For ANSWERS this effect might be attributed to the capacity based approach for water flow, because the estimation of the drainage is not completely consistent with the lower boundary condition of the lysimeter experiment.

In contrast to the water flow the modelling of the fate of MBT is rather problematic.

MARTHE, ANSWERS and MACRO account for a soil temperature and moisture dependency of degradation. The model variables soil temperature and soil moisture near the surface are well described (Table 5), whereas the degradation of MBT is not that well described. It seems like the modelling of soil moisture and temperature is easier than the quantification of the relation between degradation and soil moisture/temperature. Basically the models can be divided into two groups. One group is TRACE and MACRO using a long half-life and predicting the long-term behaviour correctly. The other group is MARTHE and ANSWERS using a clearly shorter half-life and predicting the amount of MBT residues of the first sampling date (252 days after application) correctly but for the last sampling (627 days after application) the residues are underestimated. Probably the first order kinetics approach applied for all models is not suitable to describe the degradation of MBT (Diekkrüger et al. 1995b). The first order kinetics approach might by appropriate to describe the short-term behaviour of the MBT degradation, but the long-term behaviour might be influenced by a very slow sorption/desorption kinetic. Then a hockey stick degradation function might be more appropriate. A validation of this statement from the data used in this study is limited, because a longer monitoring with more sampling dates would be required. This would also be necessary to decide if the moisture and temperature dependent first order decay is superior to the simple first order decay. The degradation rates of MBT in literature vary between 30 days [START_REF] Rouchaud | Metabolism of methabenzthiazuron in the soil of pea crops[END_REF], 139 days [START_REF] Jarvis | Simulation of soil water dynamics and herbicide persistence in a silt loam soil using the MACRO model[END_REF] and 172 days [START_REF] Wüstemeyer | Abbau und Transport von Ethidimuron und Methabenzthiazuron in einer Parabraunerde-Beitrag zur Validierung des Lysimeterkonzepts[END_REF] under field conditions. This might be another hint to a slow sorption/desorption kinetic for the longterm behaviour of MBT.

A leaching of 0.0059 % of the applied mass of MBT occurs. This process can be described with MACRO only after calibration. Without the calibration of the macropore soil hydraulic properties and the fraction of sorption sites in the macropores no leaching was estimated with MACRO. This contributes to the findings of [START_REF] Bergström | Evaluation and comparison of pesticide leaching models for registration purposes[END_REF], where only after calibration the pesticide leaching due to preferential flow is well described. For registration purposes in Europe a leaching level of 0.1% or more of the applied dose is relevant (Vanclooster et al., 2000). From this point of view the small amount of MBT leaching is not relevant, nevertheless there is a small lack of transport process description detected for TRACE, ANSWERS and MARTHE.

For the reproduction of evapotranspiration, soil moisture and drainage basically all model concepts are suitable. For the description of degradation the performance of the simple first order kinetics is comparable to the temperature and soil moisture dependent first order concepts. Here a more extensive data set on MBT residues will probably reveal the better performance of the more sophisticated approach taking soil moisture and temperature into account. Further improvement of the description of degradation might be possible with a nested first order approach allowing for different half-lifes for the short-term and the longterm behaviour of pesticide degradation. The biggest conceptual constraint concerning pesticide fate was detected for ANSWERS. Due to the thickness of the uppermost element, necessary for the Green and Ampt infiltration, the transport of the pesticide cannot be described if the downward movement of the compound is very small. The center of mass predicted with ANSWERS is always much too deep in the profile. And degradation and sorption occur usually very close to the surface, which cannot be described properly with ANSWERS. The other conceptual constraint can be detected for MARTHE, TRACE and ANSWERS. These models are unable to describe preferential flow, which must be seen against the background of the difficult parameter identification. The problem of the derivation of macropore parameters will increase with scale and a process can only be taken into account if also the parameters are available.

The use of five parallel lysimeters revealed a high variability of measured drainage amounts, while the evapotranspiration is rather similar. The cores were sampled quite close to each other (~ 1m), indicating that the spatial variability of soil properties is high. Even the use of five parallel lysimeters involves still a large uncertainty concerning the field scale drainage amount. Thus the lysimeter concept is probably not able to capture the field scale heterogeneity of processes strongly influenced by the spatial variability of soil properties.

CONCLUSIONS

During the calibration of plant parameters three crucial parameters were identified: the LAI, the root depth and the crop conversion factor. Among these parameters the crop conversion factor values were changed significantly during calibration. For most cases the K c was increased which is seen as a result of the lysimeter set up for the measurement of the actual evapotranspiration. Generally the selected databases of plant parameters are applicable to the 'Zwischenscholle' test area.

Having in mind that only the plant parameter values were calibrated, the water flow is well described with all the model concepts, which is proved by the calculated validation criteria.

The fate of MBT is not as well described as the water flow. The residues of MBT are estimated for the right depth with MACRO, MARTHE and TRACE. But the amounts of MBT residues are only poorly reproduced. Laboratory measured soil moisture and temperature dependent half-life does not improve the modeling of degradation significantly. More extensive data on MBT residues would be needed to verify a slow sorption/desorption kinetic, which also calls the first order approach applied for all model used in this study into question.

A very small amount of MBT leaching can be described with MACRO, but only after extensive calibration.

Basically all the model concepts are well applicable for water flow. The amount of preferential flow not considered with TRACE, ANWERS and MARTHE is not significant for the water flow but for the MBT leaching. A conceptual limitation was found for the capacity based ANSWERS. Here the Green and Ampt equation for infiltration requires a very thick uppermost element hindering the sound estimation of transport and degradation of MBT close to the surface. according to Graham-Bryce [START_REF] Wüstemeyer | Abbau und Transport von Ethidimuron und Methabenzthiazuron in einer Parabraunerde-Beitrag zur Validierung des Lysimeterkonzepts[END_REF] and according to Boesten and van der Linden (1991), with reference DT 50 =162 d, θ ref =0.39 cm 3 cm -3 , α= 0.01 K -1 , β=0.2 and T ref =20°C (Eq. 6). Volumetric water content was transformed to gravimetric with the topsoil bulk density of 1.57 g cm -2 .

FIGURE CAPTIONS

  account for temperature or soil moisture dependent degradation. For this case a first order kinetic decay with a DT 50 of 200d was found after calibration. By using the mean of the measured soil temperature at 10 cm depth (10T =10.3°C) and the mean of the measured soil moisture at 25 cm depth ( 25  =0.22) a DT 50 for reference conditions (T ref =20°C

Figure 4 .

 4 Figure 4. For the depth of 25 cm the first significant deviations occur during the drying period

Fig. 1 .

 1 Fig. 1. Measured atmospheric conditions and calculated reference evapotranspiration during

Fig. 2 .

 2 Fig. 2. Isoline plot of soil temperature and moisture dependent DT 50 for MBT calculated

Fig. 3 .

 3 Fig. 3. Measured and predicted actual evapotranspiration (cumulative). Bars are standard

Fig. 4 .Fig. 5 .Fig. 6 .

 456 Fig. 4. Measured and predicted soil moisture at two depths. Bars are standard deviations of

Fig. 7 .

 7 Fig. 7. Measured and predicted MBT residues (applied MBT = 100%) 252 days (a) and 627

Fig. 8 .

 8 Fig. 8. Accumulated mass of measured MBT in drainage water and model results

  a model with CME lower than 0.5 should be rejected. With a CME of 0.1 ANSWERS is below this threshold and even TRACE is quite close to this value having a CME of 0.58. The use of a criterion like the CME is supposed to support the validation of models by bringing in an objective component, but it is still the model user to decide whether a model should be rejected or not. The IA for the drainage estimated with ANSWERS and TRACE, 0.67 and 0.89 respectively, is in an acceptable range. The total amount of drainage estimated by ANSWERS and TRACE is close to the measurements. And Figure5also reveals, that the variation of drainage with time is generally reproduced. Thus the modelling of drainage with ANSWERS and TRACE is still acceptable, although this variable is not as well reproduced as the other water flow variables. It would be expected that the drainage flow could be well described, if also the soil moisture close to the lower boundary of the lysimeter is well described. This is true for MARTHE and MACRO, but not for ANSWERS and TRACE.

ACKNOWLEDGEMENT

The authors acknowledge financial support by the EU within the framework of PEGASE, contract EVK1-CT1999-00028.


Thiéry

Dominique Gutierrez

From modelling soil columns to large scale aquifers: an illustration of the MARTHE code capabilities

In the last 20 years, numerical modelling of Groundwater flow and transport has become the best tool for water resources management and decision making. Numerous codes deal with flow, transport, density, unsaturated zone, surface water-groundwater interaction, etc…but only a few codes integrate these aspect simultaneously.

The MARTHE code presented here tackles all of them in a dynamic way. This presentation intends to illustrate through case studies, how this modelling software can be used to treat various problems.

.

Introduction

Advances in computer resources observed in recent years have favoured and enlarged the use of modelling softwares for the treatment of hydrogeological problems. Memory, computer speed and storage become less and less a limiting factor. Technological improvement tend to popularize the use of numerical models. However, modelling has increased in complexity, dealing not only with flow and transport but also with density, variably saturated medium, reactive transport, river water / groundwater interactions, multiphase flow and so on.

Hydrodynamic and hydrodispersive modelling became during 20 last years, the principal management tool of the water resources. The possibility of taking into account space and temporal heterogeneities, by discretizing within a grid study area and by cutting out time periods called "time steps", makes it possible to the hydrogeologists to bring a quantitative response to problems which would otherwise only get approximate solutions. Without giving to model results an exaggerated importance, knowing of course that the value of the model depends primarily on that of the input data and of their analysis, the tool which we describe allows: to compute water balance over aquifer systems as well as specific zones (such as sub-catchment areas) defined by the user, distinguishing the amount of water corresponding to recharge, boundary conditions, storage, leakage, connexions between aquifers, and so on) to take into account the impact on the water table of civil works or hydraulics engineering (underground car park, excavations, tunnels, gravel pits, dams, etc.) to consider the interference between wells, to manage the abstraction, to work out scenarios of exploitation.

to simulate the propagation of a pollutant in a variably saturated medium, considering the characteristics of the pollutant such as degradation, sorption or density.

to represent the spanning extension of a contamination in order to evaluate its impact, the protection measures and the remediation scenarios.

to envisage the impact on the groundwater quality of specific installations such as waste disposal or industrial settlement.

For more than 20 years BRGM has worked in the field of groundwater modelling and has developed a code, MARTHE [START_REF] Thiéry | Modelling contaminant transport through the unsaturated zone in transient state with a random walk particles scheme[END_REF](Thiéry , 1993a(Thiéry , 1993b[START_REF] Thiéry | Calibration of groundwater models by optimization of parameters in homogeneous geological zones in "Stochastic and Statistical methods in hydrology and environmental engineering[END_REF] able to tackle these problems in an integrated and dynamic way. The following description of the MARTHE software will not detail the code and numerical implementations. Its purpose is to share a user experience and show the ability of the software to answer a wide range of hydrogeological problems.

Main features

The field to be modelled is discretized according to a rectangular grid of "Scottish" type: each line and each column of cells have a constant width but the user is free to adapt this width according to local heterogeneities, to the density of information and to the required precision. It is moreover possible to locally refine the grid using nested (multi-stage) grids (cf. Figure 1). These features bring to the MARTHE software a flexibility generally not accessible to other finite differences models. Input and output data files are text files which can be edited and manipulated with MARTHE interface, or with another tool (for example Excel to visualize simulated series)

Calculation of flows (hydrodynamic)

Hydrodynamic calculations are carried out according to a method with finite volumes. They make it possible to answer, for example, with the following situations:

-One-layered or multi-layered aquifers (stacking of aquifers separated by more or less continuous semi-pervious strata) or in full 3D.

-Free, confined or semi-confined water table in steady state or transient state. (cf. Figure 2) -Taking into account of discontinuities such as free surface water (lakes, gravel pits), local dewatering, (and rewetting), including in multi-layer aquifers, water table overflows (river, sources, drains), tight limits (sheeting piles…).

-Horizontal and vertical permeability anisotropies.

The resolution algorithms use methods of conjugated gradients with possibility to select the more adapted solver to the problem (Choleski, Einsenstat, Orthomin) Modelling of the lower Eocene water table between the Marne river and the Oise river around Paris, using a one-layered model. It is a simple example where the permeability of the aquifer is homogeneous and where a fine knowledge of local geology (top and bottom of the aquifer layer) allows a good restitution of the heads by the model.

Observed head distribution 1972

Simulated head

The model can represent a 2D map, a cross section, a portion in cylindrical coordinates or a 3D system. The case study #1 presents the modelling of a huge aquifer system using a 3D grid. The layers can be conformable or unconformable the ones over the others, creating short-circuits between aquiferous levels. This representation allows a perfect analogy between the layers of the model and the geological layers.

Case study # 1 : Transient modelling of the SAQ and overlaying aquifer systems (Saudi Arabia)

Km Km Km Km Km Km Km Km Km 0 0 0 0 0 0 0 0 0 500 500 500 500 500 500 500 500 500 The Saq aquifer is a huge sandstone reservoir, over 375 000 sq Km (2/3 of France total area), in the North of Saudi Arabia. This area counts about 16 000 wells out of which only 600 are used for drinking water supply, the rest being used for irrigation in the desert.

Saudi Arabia

Mass transport calculation

MARTHE is frequently used to simulate the dispersive migration of an effluent in the aquifer.

According to the type of problem arising (predominance of the convection or dispersion), the user can choose between four techniques of calculation to optimize the transport scheme implemented:

-Finite volumes; -Method of characteristics (MOC) using particles; -TVD method (Total Variation Diminishing) in transient state;

-Method of Random Walk, known also as particle dispersion method; All these methods perfectly respect the mass conservation.

Transport in aquifer is simulated under its advective, diffusive and dispersive components. Four options are proposed for reactive transport:

-Exponential decay of the effluent with time; -Filiation of sub products (metabolites); -Retardation factor using the partition coefficient kd (phenomena of adsorption-desorption).

-Freundlich and Langmuir isotherms

In the case study #2, a 3D model is used to simulate the propagation of chlorides in the alluvial aquifer of the Alsace Valley in the East of France. This modelling combines transport and density effects.

River -aquifer coupling

An hydrographic network can be integrated into the model. In this case, the bed of the river is cut out in sections of river having particular characteristics such as the permeability and the thickness of the banks, the dimension of the bottom, the width and the length of the section. The sections are superimposed on the grid. They receive water from upstream and convey it downstream. They can drain or feed the water table according to the height of the water calculated in the river and the head of the water table. Sections can be dried up temporarily.

Withdrawals or direct injections in the river can be taken into account. The coupling between the height of water in the river and the discharge of river is carried out in a dynamic way by a Manning-Strickler relation. The coupling between the river and the water table is also dynamic and can be coupled with a hydro climatic water balance assessment over a specific zone (such as a catchment area), using rainfall and potential evapotranspiration records.

Case study #3 illustrates a simulation with river -aquifer interactions. For this modelling, the river valley is defined using a nested grid with fine cells (as in Figure 1). The calibration was performed in transient state from

1989 to 2003 and validated later using updated observations. This type of model is useful for predicting the location of flood events and helps the authorities in taking proper decisions.

Calculation in the unsaturated zone

Sometimes the modelling of a natural or experimental system requires to take into account the unsaturated zone between ground surface and the water table. This taking into account utilizes characteristic laws of the terrain (retention law and permeability law) which represent the variation of the capillary pressure and the permeability according to saturation of the medium. MARTHE proposes the choice between several usual schemes. For example, with regard to the retention law, the models of Van Genuchten, Brooks-and-Corey, homographic or logarithmic curves are implemented. For the law of relative permeability, the schemes of Van Genuchten, Brooks-and-Corey, a power or an exponential law are available. The treatment is continuous between the saturated zone and the unsaturated zone.

To meet needs related to modelling for the transfer of the pesticides, the taking into account of the vegetation was integrated (cf. Figure 3). The incidence of the development of the plants (leaf and root) on the water content of the terrain, and consequently on the propagation in depth of a contaminant can be thus simulated. Research are still on-going in this field (cf. case study #4). One of the most recent developments of the MARTHE model was the integration of the macroporosity in the ground. The capillary pressure, which governs the flows in the unsaturated zone, is represented here according to time (in days) at 50 cm of depth in the ground of an experimental field. Strong suctions are created (and are simulated) when the growth of the plant is maximal and when precipitations are scarce.

Suction

Automatic calibration, optimization and sensitivity analysis

An option of MARTHE makes it possible to automatically obtain a field of permeability starting from a given piezometry, by inverse modelling. This functionality may be practical to rough a problem and to identify zones of permeability in the absence of data. The recommended methodology remains however to define zones of permeability according to geological and hydrogeologic criteria. The calibration can then be optimized by entrusting to the model the choice of the best parameters between intervals of values fixed by the user.

The results and forecasts given by the model correspond to a single response in reaction to a set of requests (inputs of the model). Various parameters of the model may come from estimates rather than measurements (for example storage coefficient or permeability which are generally only known at discrete locations). It is thus useful, even essential, to accompany any modelling by an a sensitivity analysis. MARTHE proposes an automatic method to carry out this sensitivity analysis and to possibly determine confidence intervals.

Case study #4 : Modelling the unsaturated zone. Work in the frame of AquaTerra project

Within the Aquaterra project, BRGM together with Avignon University (France), DTU, Liège and Gembloux Universities (Belgium), is studying a small catchment located 70 km North west of Paris. At Montreuil sur Epte, a spring was used for drinking water until 2001 when the authorities prohibited this use because of its level in pesticides (Atrazine 0.2; DEA 0.7 µg/l) and nitrates (65 mg/l). The catchment is studied since 2000: farmers were convinced to replace Atrazine with Acetochlor. We are following Acetochlor in the soil and waiting for its appearance at 19 piezometers and at the spring.

Works regarding modelling have been done in a progressive manner: starting from 1D column with the ambitious objective to finish with 3D model. The soil zone (1st meter) is known by numerous measurements (16 cores/plot at 5 depth and 6 different dates during year 2000) regarding water content, Bromide and Acetochlor, which were both spread simultaneously in April 2000. We consider a column constituted of 2 horizons of soils (H1 and H2) overlaying the Lutetian limestone (over 30 m of unsaturated zone). Soil properties (retention curves and permeability laws) comes from a fitting of Van Genuchten law to observed soil characteristics (using the RETC program).

Simulation of water content with depth, under several hypothesis of vegetation effect are presented here. Dots represent average of water content (16 samples) with standard deviations. Date t=0 corresponds to 01/01/2000 and observed rainfall and PET are applied as surface condition with a daily timestep. MARTHE offers the possibility to use several methods which are used in well known soil models, to take into account the vegetation (i.e. to consider the development of roots and leafs of the plant: here maize). The MACRO scheme is the closest from observed values. If we do not consider vegetation, a discrepancy appears with depth (black line). Bromide can be considered perfect tracer and proper simulation of flow is the condition of success. Although extremely simplified, the 1D modelling globally respects the concentrations and delay of Bromide transfer, 7 days after application.

When it comes to pesticide transport, new parameters have to be introduced: mainly sorption which vary with depth (linked to Carbone organic content) and degradation which vary with temperature and water content. Sorption and degradation were estimated from soil samples in laboratory, giving the above results.

It would be possible to get a better fit by playing with parameters, however the purpose was here to see whether pesticide transport might be simulated from parameters obtained by measurements (without calibration), considering the simplifications of the modelling .

Other features

The correct resolution of calculations is checked by criteria of convergence. Convergence can be controlled following the variation of average (and maximum) head between two successive iterations, as well as the residual error of flow computation. For the latter, a total residual discharge is given for each iteration, but the cell recording the strongest residual flow is indicated. It is also possible to follow convergence in one or more particular cells of the model.

Numerous options, which are not described here, are available in MARTHE:

-Calculations of fresh Water/salt Water interfaces in monolayer or multi-layer systems.

-Density effects induced by heterogeneous salinities and/or temperatures.

-Calculation of the temperature field, for geothermic problems or energy storage.

-Taking into account of the viscosity variations with temperature.

-Coupling with a drain network (for example to simulate a gallery or a karstic system).

-Coupling with water balance modelling over hydro climatic zones.

-Taking into account of vertical fractures using equivalent transmissivity.

-Transport with interactions between water, chemical species and porous matrix.

-Multiphase flow: water and gaz ; water and oil.

-Gaz flow.

The MARTHE software consists of a computing core written in FORTRAN 95, and of an interface called WinMarthe (written in Visual C++) which plays at the same time the role of pre-processor for data preparation and of post-processor for results visualization. The WinMarthe interface allows in particular to work out the grids, to display and modify the values on a cell basis, on a selection of cells or on a whole layer, to visualize the data in a given layer and in cross sections, to visualize the model using a 3D viewer, to interpolate the results, to draw path lines, to represent velocity vectors… WinMarthe also offers numerous ways of importing and exporting data to other software (such as GIS software).

Conclusion

Illustrated by a series of examples, the use of the MARTHE software, developed at BRGM, has been described. It allows a large variety of simulations, from single multilayered columns to huge aquifer systems. From a user point of view, the software and its interface, benefiting from years of improvements, can be considered as one of the most powerful tool dedicated to groundwater modelling.

Figure 1 :

 1 Figure 1: Example of a nested grid The Somme river basin: The river bed is discretized by fine cells to better represent topography and better simulate the risks of flood event.

Figure 2 :

 2 Figure 2: Observed and simulated head distribution in a monolayer steady state modelling. Modelling of the lower Eocene water table between the Marne river and the Oise river around Paris, using a one-layered model. It is a simple example where the permeability of the aquifer is homogeneous and where a fine knowledge of local geology (top and bottom of the aquifer layer) allows a good restitution of the heads by the model.

  was later simplified to 13 significant formations (7 aquifers and 6 aquitards). The spatial discretisation consists in a regular grid of square cells (2x2 km). The grid covers approximately 1000 x 1000 km. It is composed of 3 millions of cells, out of which 1.5 million are active cells. The use of satellite image of various dates allowed the reconstitution of irrigated areas and associated abstraction in the last 20 years. The transient state was calibrated with time series from about 50 observation wells. It is run with a yearly time step over 45 years.

  (NaCl) issued from exploitation, after extraction of potassium, where disposed in heaps ("terril"). Precipitation, which have leached these heaps, have contaminated the groundwater. By runoff on the banks or infiltration through the heap, meteoric water becomes brackish (350 mg/l) before reaching surface water (Rhine) or groundwater where it is diluted. Layer 1 Contamination was maximum in 1972 when 170 km 2 of the water table had a chloride (Cl) concentration over 200 mg/l. Since 1976, direct input of Cl in river have stopped. Other measured include drilling of wells to fix and treat the pollution (70 depollution wells in 2000), accelerate the salt dissolution by aspersion (photo) and recuperation of water in trenches and wells, reduce infiltration on heaps with a tight cover. Layer 2 The model was used to simulate propagation of salts under given scenarios. The above simulation results show impact of density in the propagation of Chloride in the groundwater.L., ELSASS.P. (2006) Modelling aquifer salinity in the Potash Basin (Alsace), in Proceedings of the AIH international workshop "Gestion des grands aquifères" -Dijon -France -30 April -2 May 2006. (in press) Case study #3 : River -Aquifer relationship : Flooding in the Somme river valley The flood in the Somme river valley during the winter of 2000-2001 highlighted the possible contribution of groundwater to flooding and led to observations that can be extrapolated to other regions. A flood phenomenon is the combination of rainfall, runoff and groundwater flow. The proportion of each component in this combination depends on the basin's morphology, hydrogeology, surface conditions (including frozen soil in winter), land-use and rainfall records. The Somme river valley (France): 3D representation of the 2001 flood in the south of Abbeville. The cells are represented by parallelepipeds The geological map is draped over. Fine overflowing cells are coloured in blue to represent flooded areas. The model has shown that the contribution of compartments of the chalk aquifer to the Somme discharge increased proportionally upstream to downstream. While the groundwater contribution to stream flow during the summer months is nearly 100%, it was estimated to have accounted for around 80% during the 2001 flood. References: Pointet.T., Amraoui.N., Mardhel.V., Negrel.P., Noyer.M.L., Pennequin.D., Thiery.D. -2005 La contribution des eaux souterraines aux inondations. Modélisation des hautes eaux de la Somme., Géosciences -La revue du BRGM pour une terre durable, N°2 -pp. 80-85 Pointet.T., Amraoui.N., Golaz.C., Mardhel.V., Negrel.P., Pennequin.D., Pinault.J.L. -2003. La contribution des eaux souterraines aux crues exceptionnelles de la Somme en 2001. Observations, hypothèses, modélisation, La Houille blanche, p. 112-122 N. Amraoui, C. Golaz, V. Mardhel, Ph. Negrel, V. Petit, J.L. Pinault, Th. Pointet 2002 -Simulation par modèle des hautes eaux de la Somme. BRGM report BRGM/RP-51827-FR.

Figure 3 :

 3 Figure 3 : Simulation of the capillary pressure in unsaturated zone under a corn fieldThe capillary pressure, which governs the flows in the unsaturated zone, is represented here according to time (in days) at 50 cm of depth in the ground of an experimental field. Strong suctions are created (and are simulated) when the growth of the plant is maximal and when precipitations are scarce.

  : THIÉRY D., GOLAZ, C., GUTIERREZ A., FIALKIEWICZ W., DARSY C., MOUVET C. & DUBUS I.G. 2004 -Refinements to the MARTHE model to enable the simulation of the fate of agricultural contaminants from the soil surface to and in groundwater. Proceedings of the COST international workshop, Saturated and unsaturated zone: integration of process knowledge into effective models, Rome, Italy, 5-7 May 2004, Eds. Aagaard P., Bedbur E., Bidoglio G., Candela L., Nuetzmann, G., Trevisan M., Vanclooster M. & Viotti P., Publ. La Gollardica Pavese, ISBN 88-7830-387-9, pp. 315-320.

Figure 4 :

 4 Figure 4 : The WinMarthe interface Many tools are available for the modeller to work out the data or to display the results.

Model area Model area Model area Model area Model area Model area Model area Model area Model area

water content (%) at t = 270 days (27/09/2000)

  Bromide was applied (as well as pesticides) on the soil at day 110(20 of April 2000). At 15 cm depth, the "tracer test" is well simulated. LeachP and MACRO give similar results. However at 35 cm, the best simulation is oddly the "no vegetation" simulation.

								[Br] (mg/kg) à 15 cm 12
								10
												Avg obs.	No vegetation
								8
	This shows the difficulty and sensitivity of such
	simulations.						6
	After 300 days, Bromide observed in the soil could come from decomposition of roots and leafs from		4
	plants which absorbed part of the Bromide (measurements in plants showed that up to 50% of Br	2
	was absorbed this way). The MARTHE code takes into	0
	account this consumption by plants but does not		100		150		200	250	300	days	350
	consider the plant decomposition after harvesting.
	Moreover, the soil is reworked during winter and mixing occurs in the 40 first cm.		0.00 10.00	12.00	14.00	16.00	18.00	20.00	22.00	24.00	26.00	28.00
								-0.10
								-0.20
								-0.30
								-0.40
							z (m)	-0.50				avg obs.
								-0.60				no vegetation
								-0.80 -0.70				LeachP Agriflux	vegetation	schemes
								-0.90				Macro
								-1.00
	28
	Water content (%)
	26	15/06/2000
	24				10/10/2000					The modelling is in transient state. Vegetation effect becomes important
										during summer time. AGRIFLUX
	22									scheme	was	probably	badly
										parameterised	because	the	plant
	20									development is too fast. LeachP and
	18					avg obs. no vegetation		MACRO schemes are giving similar results. If no vegetation is considered, the
	16					Agriflux Macro			water content during summer time is over estimated.
	14					LeachP		days
	100	150	200	250	300	350	400	450	500

Vol.