Estimating the degree of preferential flow to drainage in an agricultural clay till field for a 10-year period

The conceptual understanding of the preferential water flow is crucial and hence understanding the degree of water percolating rapidly through vertical macropores, or slowly through the low-permeable matrix, is vital in order to assess the risk of contaminants like nitrate and pesticides being transported through a variably-saturated macroporous clay till to drainage. This study compared six different model concepts, using the dual-permeability module of the one-dimensional model DAISY, incorporating three different macropore settings and two different groundwater tables set as lower boundary conditions. The 5 three macropore settings included vertical macropores supplying water directly to a) drainage, b) drainage and matrix and c) drainage and matrix including fractures supplying water to the matrix in the saturated zone. The model study was based on ten years of coherent climate, drainage, and groundwater data from an agricultural clay till field. The estimated drainage obtained with the six model concepts was compared to the measured drainage. No significant discrepancies between the estimated and measured drainage were identified. The model concept with the macropore setting b) exposed to groundwater fluctuations 10 measured in the southern part of the field, gave the best description of the drainage. Bromide leaching tests were used to evaluate the mass balance of the model concepts. The estimated water balance of all six concepts revealed that 70% of the precipitation input to drainage was transported via macropores. According to the results of bromide leaching simulation, 54% of the drainage was estimated to be transported via vertical macropores being initiated in the plow layer.


Introduction
Preferential flow (PF) is a general term for non-uniform, non-equilibrium water flow processes, which can occur in low permeable drained clay tills given the presence of macropores/discontinuities/continuous voids such as earthworm, root channels, sand lenses, and fractures. Today, it is well-recognized that such macropores can provide the main transport pathways for leaching of agrochemicals and nutrients from the soil surface to drainage and groundwater. The occurrence of PF and the 20 associated solute transport have been observed and demonstrated under several field and laboratory conditions (Akay and Fox, 2007;Gjettermann et al., 2004;Larsbo et al., 2009;Nielsen et al., 2010;Petersen et al., 2012;Rosenbom et al., 2009;Stamm et al., 2002;Tofteng et al., 2002). Studies have documented that drainage is affected by the presence of macropores, especially 1 earthworms (Akay et al., 2008;Frey et al., 2016;Kohne et al., 2006;Larsson and Jarvis, 1999;Petersen et al., 2012;Vogel et al., 2000) providing a direct hydraulic connection between the soil surface and the tile drain. This direct connection may 25 result in rapid transport of nutrients (Cheng et al., 2014;Gardenas et al., 2006;Kohne et al., 2006;Larsbo et al., 2009;Larsson and Jarvis, 1999;Mohanty et al., 1998) and other contaminants from the surface into the tile drain system. Bypassing the soil matrix and its retarding capacity (sorption and degradation), the contaminants are transported via tile drains to nearby surface water bodies, such as lakes and rivers, with no or limited possibility of further reduction or retention along this path (Blann et al., 2009;Hansen et al., 2012a;Hinsby et al., 2012;Kronvang et al., 2005). This direct connectivity phenomenon was 30 verified by Akay and Fox (2007), who conducted infiltration experiments in soil columns with an artificial macropore placed directly above or shifted from the subsurface tile drain pipe. The set-up of their experiment allowed the length of the artificial macropores to be varied without disturbing the soil column. The macropores were classified as either surface-connected or buried macropores. They observed that the water arrival time to the tile drain through surface-connected macropores was eight times faster, and about two times faster for buried macropores, than in a soil column with only matrix flow. At field scale, Pe-35 tersen et al. (2012) showed a direct connection of macropores to the subsurface tile drains by injecting smoke to the tile drain pipes and counting the surface-emitting macropores on sandy loam. They analyzed two 35-m-long tile drain sections for the spatial distribution of emitting macropores while the water table dropped from about 0 cm to 20 to 50 cm below drain depth.
From 65 to 96 smoke emitting macropores were identified along the tile drain sections, separately. Geologies incorporating such highways for flow and transport in the shape of macropores connecting the surface/plow layer with the deeper part of the 40 soil can, as described, make them vulnerable in regard to rapid leaching of contaminants to drainage and groundwater. It is, therefore, imperative to be able to account for these highways when assessing the risk of leaching of contaminants through the variably-saturated drained clay till. As demonstrated by Germann (2018) this is not straightforward, since the hydromechanics of the macropores and the interaction with the lower-permeable soil matrix still pose a scientific challenge in terms of avoiding too crude a description thereof. Nevertheless, within the last four decades, different numerical models have been developed, 45 incorporating reality-simplified processes/concepts describing PF and coherent solute transport. With the macropore setting being very complex and varying in space and time, the well-documented models like RZWQM2, MACRO, HYDRUS1D and DAISY (Hansen et al., 2012b;Larsbo and Jarvis, 2003;Ma et al., 2012;Simunek and van Genuchten, 2008) cope with the water flow and solute transport in one dimension using the Richards (1931) equation for capillary flow in the matrix domain and apply different approaches for the macropore domain. Hence, in these models, the soil pore space is divided into two 50 domains: A low-permeable matrix domain with a relatively large storage capacity and a high-permeable macropore domain with a relatively low storage capacity. The Root Zone Water Quality Model, RZWQM2 describes the macropore flow capacity using Poiseuille's Law based on gravity flow in cylindrical pores. The water exchange between the domains is based on the radial or lateral Green-Ampt equation, that accounts for the infiltration front in the matrix surrounding the cylindrical pores or cracks, respectively. Using pore-scale physical laws in RZWQM2 poses a disadvantage, in that it assumes full saturation in 55 macropores and simplifies the geometry of macropores, which is usually tortuous (Simunek et al., 2003). The governing equation for the macropore domain of the dual-permeability model in MACRO (Larsbo and Jarvis, 2003) follows the kinematic wave approach described by Germann (1985), where the macropore flow is controlled by a kinematic exponent, which reflects the macropore size, distribution, tortuosity and ignores capillary flow. In MACRO, the mass transfer between domains in the absence of gravity is described by the Richards equation, recast as a diffusion equation for a homogeneous medium, where the 60 water content gradient is the driving force. In HYDRUS 1D (Simunek and van Genuchten, 2008) uses the approach of Gerke and vanGenuchten (1996) who applied the Richards equation for both matrix and macropore domain. This dual permeability approach requires more, rarely available, specific knowledge about water retention and hydraulic conductivity functions for both matrix and macropore domains and their interface. The mass transfer between the two domains is similar to MACRO, using the first-order mass transfer coefficient and assuming it is proportional to the difference of the pressure head between the 65 two domains. Also, HYDRUS 1D has incorporated different approaches to simulate the flow in the macropore domain, such as the mobile-immobile and dual porosity model, which requires fewer input data. DAISY has been modified by Mollerup (2010) based on the laboratory experiments of Tofteng et al. (2002) and Gjettermann et al. (2004) and the field experiment of Petersen et al. (2012). Tofteng et al. (2002) showed, with an artificial macropore attached to the bottom of a sandbox, that macropores were able to conduct water on pressure potentials as low as -30 cm, which indicates that macropore flow can alter from a typi-70 cal retention/conductivity function. These findings are essential, since tillage increases the porosity of the plow layer (Petersen et al., 1997), while macropore continuity may be reduced (Roseberg and Mccoy, 1992), and a compacted layer -the plow pan, which often occurs at the transition zone between the plow layer and the subsoil -may cause accumulation of water and thus initiate macropore flow (Frey et al., 2012). All of the models mentioned are simplifications of reality and as such, have pros and cons in regard to their description of the hydromechanics of the macropores and the interaction with the surrounding matrix 75 and drainage. Additionally, they need to be able to account for plant cover and fate processes (mineralization-immobilizationturnover, sorption, degradation, effects of soil management, etc.) in an assessment of the leaching risk of agrochemicals. The purpose of this study was to evaluate whether it is possible with the DAISY model to capture the water balance (climate, groundwater fluctuations, water content in the profile and in drainage) of a well-described agricultural clay till field that is tile-drained and prone to preferential transport (Ernstsen et al., 2015), during a 10-year period with/without macropores and, if 80 so, to estimate the degree of macropore flow and the associated transport of bromide to the drains during the period.

Field characteristics
In this study, data originated from a tile-drained agricultural clay till field included in the monitoring program PLAP (Pesticide Leaching Assessment Programme of Denmark) (Lindhardt et al., 2001). The field is situated at Silstrup in the North Western 85 part of Denmark (56 • 55'56"N 8 • 38'44"E). The field was established in 1999 and has been used for agricultural purposes since 1942 and systematically drained since 1966. The field has been farmed from 1983 without any plot experiments conducted (Lindhardt et al., 2001). The field and the surrounding area are dominated by Weichselian glacial clay till deposits in combination with thin meltwater sand lenses. The underlying geology comprises dislocated slices of clay till and meltwater sand/gravel forming "drains" down to the groundwater. The southern rim of the field contains dislocated clay and silt mixtures in the depth 90 interval of 4-8 m. This material is underlying half of the field (Fig. S1). According to the geological survey of the field, the top layer of this mixture is the most permeable part of the whole field and seems to dominate the hydraulics of the field-setting (Abiltrup, 1999). In general, the field is heavily fractured and contains a high-density of macropores. In the upper 2 meters of soil, dominated by bioturbation and desiccation cracks, the macropore density was estimated to be 400 pores m −2 (Lindhardt et al., 2001) with diameters of 2-5 mm, using the methodology from Klint and Gravesen (1999). The size of the field is 1.69 95 ha (185 × 91 m). Highest points of the field vary from 41 to 45 m above sea level (m.a.s.l.) (Fig. S2). The groundwater head potentials (in m.a.s.l.), measured in February 2000 from several piezometers and monitoring wells surrounding the field (Fig.   S3), indicated that the groundwater tends to flow towards the South West, which is opposite to the topographical slope of the field. This water movement suggests that the significant part of the field is hydraulically connected to the layer with the silty-clay Oligocene deposit, which can be considered as an aquifer, overlaid by a clay till (Freeze and Cherry, 1979) (Fig. S1).

100
The pedological profile investigation was carried out in the North Eastern corner of the field where two profiles were excavated (N-North face, E-East face). The most common horizon subdivision was A p -B v -B t(g) and A p -BC-C c , with an average depth of A p (plow layer) and B v /BC of 31 and 140 cm, respectively. The textural analysis showed a plow layer with 18-25% clay and more variable clay content in the underlying soil of up to 43% clay in the eastern profile. This clay content probably decreases towards the southern part of the field since Hansen (1976) found 15% and 24% clay in the topsoil and the underlying soil, 105 respectively. The field was also mapped with a DUALEM21 ground conductivity meter (GCM). The processing and inversion of GCM data were done by Aarhus Workbench using the Aarhus Inv inversion code (Auken et al., 2018) (Fig. S4). The soil hydraulic measurements for soil water retention (SWR) and hydraulic conductivity were conducted on 100-cm 3 and 6280-cm 3 soil cores, respectively. The cores were sampled at three depths corresponding to the A p -, B-, and C-horizon. The mean retention points showed that the water content gradually decreased with pF (logarithm of the soil water potential [ψ, in cm H2O 110 times -1]) because of the relatively uniform pore size distribution, which is consistent with loamy-clay soil types (Table S1).
The SWR, saturated K s and unsaturated K(ψ) hydraulic conductivity were determined by the methods described in Iversen et al. (2004) and Iversen et al. (2011). The unsaturated hydraulic conductivity showed high variability at -10 cm H2O pressure potential and was significantly lower than the saturated conductivity. The variability in saturated hydraulic conductivity showed differences by two to three orders of magnitude, which indicates the presence of macropores (Fig. S5) 115 2.2 Field-scale monitoring and field management The tile drain system of the field consists of five parallel tile drains at approximately 110 cm depth, 18 m apart, established in a south to north direction. They are connected with a transverse tile drain at the northern end of the field. Drainage flow in the drainage outlet well and soil water content were measured hourly and automatically from 2000 to 2010. The drainage flow was calculated from a measured water level above the Thomson weir V-notch (Lindhardt et al., 2001) (Allen et al., 1998). The wind speed measurements at 10 m height were scaled to 2 m, assuming logarithmic wind profile and neutral atmosphere (Allen et al., 1998). Hourly precipitation [mm] was sampled with a tipping bucket gauge (Lambrecht meteo GmbH, Göttingen,Germany) at 125 the upper North Eastern end of the field (Fig. S2).
[ Figure 1 is about here] In Figure 1, the hourly recorded drainage presented in a seasonal accumulation, where a season comprises the period from 1 st April until the following 31 st March, in order to cover an average crop growing period. The seasonal precipitation sum varied from 691 mm to 1146 mm, with an average of 910 mm for the above-considered time interval. Most of the precipitation 130 records (96%) were obtained from the field, but where the gauge malfunctioned (4% of the records) precipitation was obtained from the nearby meteorological station. The precipitation was obtained in open field conditions, 1.5 m above the soil surface and was corrected to true ground based on factors provided by Allerup and Madsen (1979). This correction was conducted in order to overcome induced errors, mainly due to wind turbulence, which caused the 1.5 m precipitation to be underestimated.
To identify the possible solute transport routes, 30 kg ha −1 Potassium Bromide (KBr) was sprayed on 22 nd May 2000. Water 135 samples were collected flow-proportionately at the drainage outlet (Ernstsen et al., 2015) for in-lab determinations of bromide (Br -) and nitrate (NO 3 ) concentrations. The accumulated bromide transport is presented for the following two seasons (season -2000, 2001) after the application in Figure 1a. Soil water content was measured with horizontally installed 30-cm-long TDR (Time Domain Reflectometry) probes at various depths (0.25, 0.6, 0.9, 1.1, 1.9 and 2.1 m with three replicates per depth) in the variably-saturated zone, where the pedological profiling had been conducted (Lindhardt et al., 2001). In this study only 140 the results of replicates at depth -25 cm and -60 cm are presented (Fig. 1b). The measurements were obtained using probes made at Aarhus University and a Campbell CR10X data logger controlling a Tektronix 1502C cable tester (Tektronix Inc., Beaverton, OR, USA). The hourly measured water content (three replicates) were averaged to daily response data, due to the hourly measurements containing a significant amount of missing data due to instrument failure. Further, as the water content measurements were sampled in the northern corner of the field, these measurements may not represent the water flow dynamics 145 and variability of the whole field and especially not the PF (Nimmo, 2012) as the 30 cm long probes integrate by far the size of a macropore, and water flow in bigger pores happens so fast that it may not be captured by hourly measurements. Germann (1985) reported velocities of the wetting shock fronts in the range of 0.3 to 1 mm s −1 . However, the hourly measured drainage flow was kept at its hourly time resolution, as it integrates the whole field response to measured hourly rainfall data. The hydraulic response data also included daily groundwater table measurements in piezometers situated at the northern (P4) and 150 southern (P3) end of the field (Fig. 1c, Fig. S2). P4 was placed in clay till in the north, whereas P3 hit dislocated clay-silt layers ( Fig. S1 and S2). With the earlier mentioned difference in geology from south to north, a contrast of up to 1.5 meter in the groundwater level was measured during the summer period and, hence, slower drainage with a lower amount may have taken place in the north compared to the south. Unfortunately, this cannot be documented by monitoring data, since only TDR has been deployed at the northern end of the field. In the course of the data period of 2000-2010, five different types of crops 155 were grown: fodder beet (2000,2008), spring barley (2001, 2005, 2009), pea (2002), maize (2003), winter wheat (2004, 2007) and winter rape (2006). All crops were harvested and analyzed for dry matter yield and nitrogen content, and the phenological development stages of the plants were recorded.

The DAISY model
DAISY is a soil-plant-atmosphere system model incorporating modules describing plant water uptake, plant growth, C-N 160 turnover, and water (drainage transport and percolation to the groundwater), heat, snow-melt and nitrogen dynamics in the soil (Hansen, 2002). Each of these processes is incorporated into modules that are interconnected with each other. In this study, the main focus will be on the soil water module (Hansen et al., 1990;Hansen, 2002;Hansen et al., 2012b;Mollerup, 2010). DAISY uses the Richards equation (Richards, 1931) to solve the one-dimensional water flow in a variably-saturated subsurface porous medium. In DAISY, the Richards equation includes the closed form of the van Genuchten SWR model 165 (Van Genuchten, 1980), which is coupled with the Mualem hydraulic conductivity theory (Mualem, 1976) applied to the matrix domain. In terms of soil solute balance, the model applies the convection-dispersion equation with the ongoing transformation processes and uptake by plants. For the modeling of PF, DAISY assumes two flow regimes: a matrix flow regime where the convection-dispersion equation for solutes is applied, and a macropore regime where only convection is considered. The mass exchange between the two regimes is governed by the water flow through the sink-source strength of the two domains and the 170 model does not consider storage of solutes in the macropores; hence the macropores are only considered as a fast pathway for transport of solutes (Hansen, 2002).The fast flow domain in DAISY is described by a macropore module described by Mollerup (2010) (web address:https://daisy.ku.dk/publications) and tested in technical reports, prepared for and published by the Danish Environmental Protection Agency (Hansen et al., 2010a(Hansen et al., , b, 2012c. Here a macropore is a vertically oriented feature, characterized by physical properties such as length, diameter and density. The macropores are divided into two classes 175 specified by their lower boundary. One class of the macropores ends in the soil matrix, while the other class can be set to end directly in the tile drain. When water is transferred from the matrix to the macroporous domain ending in the matrix, the water is instantaneously moved to the top of the current water level in the macropores. If the pressure difference between the macropore and the matrix exceeds a predefined soil water pressure barrier potential (ψ barrier ), water will be transferred back to the matrix from the macropore. For the macropores ending in the tile drain, the lower boundary condition has a constant 180 pressure head of zero, meaning that water will be transferred directly into drainage without building up in the macropores.
The macropore flow is initiated when the matrix potential exceeds a specific pressure potential called ψ init . Hence, if this pressure potential is exceeded, the macropore domain is activated and water starts to fill up the macropore. When the pressure potential drops below a level called ψ term , inflow to the macropore is terminated. All pressure threshold parameters apply to all macropore classes. The macropore domain concept consists of an equidistant distribution of macropores where the water 185 flow in the domain is quantified by soil water extraction from the matrix, which is based on the water movement in a confined aquifer, towards a well (where the well is represented by a single macropore). The interaction between macropores and surface water follows Poiseuille's Law (Hillel, 1998) in the same way as in the RZWQM2 model (Ma et al., 2012). The only driver of the water flow is gravity. The macropore model does express viscous flow, as Germann (2018) suggested, since the water flow is lumped to an instantaneous transport within a macropore. Due to the velocities of wetting shock fronts, which is the range of 190 0.3 to 1 mm s −1 , the model could not represent the hydro-mechanical aspects of a turbulent flow, since the basic time step is one hour. When the hydrologic set-up includes a subsurface tile drain at a certain depth, drainage from the soil matrix domain will only be initiated when the groundwater level rises above this depth. This inflow to the tile drain assumed installed in a flat landscape, is described by the Hooghoudt equation (Hooghoudt, 1940).

Model concepts 195
The soil profile in the one-dimensional DAISY model was described based on the pedological survey of the Silstrup field with three different horizons: -A-horizon -Plow layer to a depth of 31 cm -B-horizon -Layer from 31 cm to 140 cm depth -C-horizon -Layer from 140 cm to 500 cm depth 200 The discretization of the profile was into 1-cm layers to 50 cm depth, 5 cm to 110 cm depth, and 10-cm intervals to 500 cm depth, including additional steps at the upper and lower boundaries of the horizons. The tile drain was set at a depth of 110 cm and 18 m apart, according to the field characteristics described above. The C-horizon was extended to 500 cm in order to account for the fluctuating groundwater table. To avoid having an unrealistic impermeable layer as a lower boundary condition as the Hooghoudt equation requires, a Dirichlet boundary condition defined by the measured groundwater level was 205 applied to the Hooghoudt equation (Fetter, 2000). The van Genuchten-Mualem (vGM) retention conductivity model was used to describe the water retention and conductivity in all three layers for the matrix domain. The retention parameters (θ s , θ r , l, n, and α) of this model and saturated hydraulic conductivity K s [LT −1 ] were estimated by fitting the vGM model to the laboratory-measured water retention. The uncertainty boundaries were established by fitting all retention and conductivity data from the corresponding samples of the horizons and combining them to create a retention and conductivity range. (Table 1 and 210 Fig. 4, 5). The unsaturated hydraulic conductivity was fitted to the measured K(ψ)available from different soil cores (Table   1, Fig. 5) for each horizon. The outcome of the fitting gave saturated hydraulic conductivities that were one to two orders of magnitude lower than conductivities measured directly on the undisturbed soil cores (Fig. S5). The initial SWR and hydraulic conductivity parameters were obtained from the best fit of the combined vGM models (Fig. 4,5). In Table 1, the columns Min. and Max. show the range of the uncertainty boundaries of the parameters assessed from the fitted results. To incorporate 215 PF, a well-connected macropore system needs to be included in the model. To represent the macropore system in the field, three different macropore settings (MSETs) have been included in the model (Fig. 2, Table 2). The division into the different MSETs was based on a field study of Nielsen et al. (2010) conducted on a similar macroporous, sandy-clay loam soil. Because the macropore distribution in a given field is unknown and difficult to describe by soil textural parameters (Lamande et al., 2011), the following conceptual considerations have been applied. Matrix macropore 1 (MM1) represents all small macropores 220 created by bioturbation or agricultural processes after plowing. This kind of macropore always ends at the bottom of the plow layer. MM2 and MM3 represent the wide (up to 5 mm in diameter) and long macropores created by large earthworms or old root channels. They are present below the depth of the tile drain facilitating, in most periods of the year, a direct connection to the groundwater table. The difference between MM2 and MM3 is the starting depth. MM2 starts at the soil surface, whereas MM3 starts at the bottom of the plow layer. This difference defines a surface-connected macropore as MM2, and a buried 225 macropore as MM3, consistent with Akay and Fox (2007) and Rosenbom et al. (2009). Matrix macropores are able to increase the fast storage capacity in a clay till which can slow down the breakthrough of the solute transport. The macropores ending in the tile drain (DM1, DM2) are numbered according to the starting depth. The field description (Lindhardt et al., 2001) indicates that the field is heavily fractured, which is the reason for including a long macropore passing through the system to represent a well-connected fracture system of a field in a one-dimensional model (Fig. 2, Table 2).

230
[ Figure 2 is about here] The concepts represent a clay till profile with: a) only drain-connected macropores (DM), b) macropores ending up in the drain (DM) and the matrix (MM), and c) DM and MM macropores where an MM-type is present through the entire soil column, representing a fracture (FR). With the matrix-macropore-concept applied in DAISY, water cannot percolate out of the bottom of the MM macropore but only through its sides, while water can build up inside the macropore if the surrounding matrix has 235 sufficiently low permeability. Given this model constraint, a thin, permeable sand layer (even though not present) needed to be added at the bottom of the model profile to facilitate constant drainage of the FR macropore. The soil texture properties of this sand layer "D-horizon" (Fig. 2) was taken from the soil database at the Department of Agroecology "Jordprofildatabasen" (Madsen et al., 1992), and the vGM parameters were generated by HYPRES (Wosten et al., 1999). The initial macropore parameters (Table 2) were extracted from the field study of Nielsen et al. (2010), where they counted the stained macropores 240 along a drain trench, after a tracer experiment. MM1 density ρ M M 1 was set to 100 m −2 to represent the large amount of bioturbation. The remaining macropore properties were set arbitrarily since they were considered to have been applied in the model calibration process. The initial fracture diameter and distribution (d F R and{ρ M M 1 ) were taken from Rosenbom et al. (2009) and Lindhardt et al. (2001), respectively. To apply realistic boundary conditions for the model, hourly precipitation, evapotranspiration and crop developments were used as the upper boundary condition, and the depth to groundwater table was 245 used as a lower boundary condition. The latter was monitored on a daily basis. Given the large difference in groundwater table level monitored in piezometer P4 in the northern part of the field and P3 in the southern part ( Fig. S1) during the summer, the flow through the profile exposed to both conditions was evaluated for all three MSETs; a), b), and c): In such a clay till it is imperative to account for the fluctuating groundwater table and coherent drainage to be able to account 255 for the water balance in the system over time addressing the field, including the degree of macropore flow. This study should be seen as a "learning" study, which can help to delineate the degree of preferential flow and transport through fields in which one can only have knowledge about the fluctuations of the groundwater table and knowledge about the presence of tile drains.

Model performance and objective function
In order to calibrate the six model concepts towards the desired performance by fitting monitoring data, an automatic calibration 260 procedure was set up. Efficiency criteria were defined as mathematical measures to evaluate model performance, e.g., how well a given model fits the available observations (Krause et al., 2005). Several studies have discussed which efficiency criteria to use in a hydrological modeling context (Krause et al., 2005;Muleta, 2011;Singh et al., 2016;Willems, 2009). According to Krause et al. (2005), most of the efficiency measures that use the squared deviation tend to be sensitive to peak flows, which is an inherent characteristic in tile-drained soils. In this study, the Normalized Mean Absolute Error (nMAE) was used 265 to avoid giving extra weight to peak drainage. This performance measure was chosen with reference to Muleta (2011), who tested efficiency measures based on the minimization of the squared, absolute, and log of the error terms. For proper automated calibration, several different objectives might need to be considered in order to avoid pushing errors into an unexamined domain of the model. (Madsen, 2000) The objective function for the calibration procedure in the present study synthetizes four different sub-objectives. Primarily the hourly drainage and cumulative drainage observations were taken into account in order to take 270 control of the dynamical behavior of the model and the water quantity leaving the clay till profile. To further constrain the model and to have a better representation of the soil water dynamics, the daily soil water content represented by TDR measurements at 25 cm (SWC25) and 60 cm depth (SWC60) were added to the objective function as two sub-objectives. Since the model was constrained to different objectives with different units, the four MAEs represented in Eq. (1) were normalized by means of the corresponding observations, in order to aggregate them into one objective function, to be used for automated calibration.
,where N is the number of observations of a given objective within a year and ,where K is the number of calibration years Given the substantial differences between years, the nMAE was calculated for 280 each objective on a yearly basis. The mean of the yearly nMAEs was used as a measure of the performance of the model during the years included in the calibration period. Thus, the final multi-objective function was expressed as: For the evaluation of the models, nMAE (range [0, +∞]) was used primarily, as it is the inverse version of Volumetric Efficiency (VE) proposed by Criss and Winston (2008 (Gupta et al., 2009) are also presented. According to Singh et al. (2005) and Hansson and Hokfelt (1975), nRMSE and nMAE values, less than half the standard deviation of the measured data may be considered low, and that either is appropriate for model evaluation. In terms of KGE, a model result higher than 0.5 is considered satisfactory.
These values further called threshold(TH). chosen to cover a broad range of yearly precipitation characteristics, with marked differences in total amount and seasonal distribution. (Fig 1a). The Morris sensitivity screening method (Campolongo et al., 2007;Morris, 1991) was applied to the model 295 parameters that were considered to have an impact on the multi-objective function. The sensitivity screening was done based on the different sub-objectives. Each parameter was screened against each sub-objective separately and selected as sensitive based on either one, two, three, or all four sub-objectives. The non-sensitive parameters were set to their nominal values. The sensitive parameters were calibrated by the DEoptim algorithm. The Differential Evolution Optimization (DEoptim) (Storn and Price, 1997) was used to calibrate the model. DEoptim is a global optimization genetic algorithm, which uses biology-inspired 300 operations of crossover, mutation, and selection on a population, in order to minimize an objective function over the course of successive generations. For the search of the global optimum, the multi-objective function was used. The sensitivity analysis and the calibration were carried out with the R statistical programming language (Team, 2017). The RDAISY toolbox (Jabloun et al., 2014), with the combination of the "sensitivity" package (Pujol et al., 2017) and "DEoptim" R packages (Mullen et al., 2009) was used to perform the sensitivity screening and the model calibration.

305
3 Results and discussion

Simple calibration of the DAISY crop parameters
The dynamics of soil water content and the matrix-macropore interface depend not only directly on soil physics, but also on evaporation and plant water use. The actual evapotranspiration was determined by the reference evaporation, the crop coefficient for the corresponding vegetation and the water transport capability of the given crop root system. Therefore, the 310 crop module of DAISY was calibrated to the measured phenological development stages and harvested dry matter yield (Fig.   S9). According to Plauborg et al. (2010), by matching the simulated dry matter yields with the measured ones, the simulated evapotranspiration may be at the right level for Danish conditions, as the standard crop parameters included in the DAISY crop database are assessed from experiments in a temperate climate and based on the strong relation between accumulated transpiration and biomass production (Hsiao et al., 2007;Steduto et al., 2009Steduto et al., , 2007. However, a check of the most essential 315 parameters for transpiration such as crop development, was carried out, as these parameters may show differences between genotypes within species. Further details on crop calibration can be found in Hansen et al. (1990), (Hansen, 2002) and (Hansen et al., 2012b).

Impact of direct drainage macropores on model performance
The performance of the model applying the groundwater table measurements from P3, before (single permeability model) and 320 after the introduction of macropores (the dual-permeability model, MSET b)) without calibration of parameters, was tested.
The two model set-ups yielded 1.09 and 1.37 nMAE, respectively, for the drainage sub-objective, when tested for the drainage year 2006-2007. Even though the single permeability model "mathematically" performed better than the dual-permeability model, it was not capable of capturing the rapid drainage dynamics (Fig. 3). This is in agreement with the findings of Akay et al. (2008).

325
[ Figure 3 is about here] The non-calibrated dual-permeability model MSET b) + P3 was parameterized with the initial parameter values described in Tables 1 and 2. The model was able to depict the drainage dynamics, but the proportion of water leaving the clay till profile via drainage was almost twice that of the measured (Fig. 3).

330
For the model set-up of MSET a), b) and c), 26, 32, and 35 parameters, respectively, were identified as probable sensitive parameters (Table 1 and 2). The difference in the number of parameters was due to the macropore parameters specified for each macropore setting. Besides the macropore parameters of the MSETs, all parameters related to the SWR and hydraulic conductivity were kept for horizons A, B and C as well as the thickness of the A-and B-horizons. The latter two were included, based on the variation in horizon thickness from the pedological profile description. In addition, Table 1 presents the uncertainty 335 boundaries (columns Min. and Max.), which were used in the sensitivity analysis. For the A-horizon, the uncertainty range was extended for the vGM parameters, after the inspection of the SWR range from PLAP (Iversen et al., 2004;Lindhardt et al., 2001). By inspecting the resistivity map of the field, the corner where the sample was taken has the highest clay content in the whole field (Fig. S4). The macropore parameter ranges were based on the macropore studies of (Nielsen et al., 2010;Rosenbom et al., 2009). According to the Morris sampling strategy, 270 to 360 parameter sets were obtained depending on 340 the MSET. As previously mentioned, in order to identify the important parameters, the sub-objectives (Drainage, Cumulative Drainage, SWC25, and SWC60) were separately screened. The results of the sensitivity screening were turned into the Morris distance( ) (Ciric et al., 2012;Jabloun, 2015), which represents the Euclidean distance of the parameter from (0,0) on the µ * − σ coordinate system (Campolongo et al., 2007). Thus, all the sensitive parameters that have an impact on either one or more sub-objectives were considered for the calibration procedure. The sensitivity criteria for Drainage, Cumulative Drainage, 345 SWC25, and SWC60 were set to ≤ 0.2, 0.2, 0.1, and 0.1, respectively. The criteria were set arbitrarily by observing the sensitivity screening outputs (Table 1 and 2, Fig. S6-S8). The sensitivity screening gave values of 14/16, 13/15, and 19/19 for model with MSET a), b) and c) with lower boundary P3/P4, respectively. All non-sensitive parameters were set to their nominal values and excluded from the calibration.
[ Table 1 is about here]

350
The sensitivity screening showed that for the models, the A-and B-horizon water retention and hydraulic conductivity parameters dominated (Table 1, Fig. S6-S8). This was anticipated, due to their strong influence on the sub-objectives (SWC25 and SWC60). In the macropore domains, the pressure parameters were also found to be sensitive in all models, with the exception of ψ barrier in model MSET b) -P3 ( Table 2). The most interesting result of the sensitivity screening of all the models was that the MM dead-end macropores had almost zero influence on the flow through the profile. This is probably because the MM 355 macropore is filled up with water quickly without being able to release water through the low permeable surrounding matrix at the bottom of the macropore. The MM macropores contribute very little thereafter to the rest of the hydraulic environment.
Neither did the diameter of the macropores show any effect on any of the objectives, which is consistent with the study of Ahuja et al. (1993), who stated that to determine the maximum flow rate, one only needs to adopt a relevant pore size, instead of identifying the distribution of the macropore sizes in a given field.

360
[ Table 2 is about here]

Calibration of the water balance of the six models
In over 30000 model runs, the multi-objective function (nMAE) was minimized to 0.75-0.85 for all six models using data from the period 1 st April 2003 to 31 st March 2008. The calibrated sub-objectives are shown in Table 3. The outcomes of the calibration of MSETs a), a), and c) did not differ significantly from each other, although, according to KGE settings, a) 365 performed significantly better than the other MSETs (Table 3). The more complex MSET c) (incorporating a fracture together with the D horizon) underperformed (nMAE=0.81-0.85) compared to MSET a) (nMAE=0.79/0.81) and a) (nMAE=0.75/0.79).
This may be because the fracture concept cannot be incorporated in DAISY; not being in direct contact with the tile drain has no significant impact on the drainage. The different lower boundary condition (P3 and P4) had little influence on the calibration outcome, but the performance measures disagreed on which of them had the greatest impact on the modeled field. If nMAE is 370 taken into account, P3 outperforms P4 in all model concepts; for KGE and nRMSE, the exact opposite is true. According to the TH, the multi-objective resulted in satisfactory results in all models according to nMAE and KGE, but non-satisfactory by nRMSE.
[ Table 3 is about here] 3.4.1 Soil water retention and hydraulic conductivity.

375
The calibrated hydraulic parameters for the A-horizon differed in all models compared with the initial SWR range from PLAP ( Fig. 4). The calibrated models showed a more realistic representation of sandy clay loam, while the SWR range from PLAP showed a very high air entry pressure, which is more characteristic for clay soil. The constant water content at the wet end was reduced to higher water pressure from -50 to around -4 cm, reflecting conceptually that air entry pressures in all models were increased, thereby increasing the potential for infiltration and redistribution within the plow layer. The SWR of the B-380 and C-horizons remained within the observed ranges. In order to validate the calibrated SWRs of the A-horizon, the observed SWR ranges from previous studies conducted on this field (Katuwal et al., 2015;Naveed et al., 2016;Norgaard et al., 2015), where 65 soil samples were taken from the whole field in a 60 x 165 m grid, are superimposed in Figure 4 (orange band) over the data from PLAP and the calibration. As Figure 4 shows, all calibrated SWRs mimic the shape of the SWR range of the studies.

385
[ Figure 4 is about here] The calibrated hydraulic conductivities were all within the observed ranges (Fig. 5), although the C-horizon conductivity was at the higher end of the observed range. For the applied two different lower boundary conditions (P3, P4), the MSETs showed similarity between the calibrated models in SWRs and the measured, unsaturated hydraulic conductivity.

Soil water content at 25 and 60 cm 390
The SWC measurements at 25 cm depth were well described by the six models ( Fig. S10 and S11, upper graph). The models captured the groundwater fluctuations and plant water uptake in the growing season, which is causing drying in the A-horizon, and the models reacted to the sudden saturation from the surface due to precipitation. There are discrepancies compared to the measured values due to the expected heterogeneity of the drained field not being taken into account by the TDR point measurements. The ability to describe the SWC measurements at 60 cm depth was not as good as for the upper depth, due to In the drainage year 2006-2007 (part of the calibration period), the dynamics of the drainage were captured completely when applying the three MSETs (Fig. 6), although they did not capture timing and some of the amplitude of the measured drainage in March 2007. In this period, the temperature was around or below zero degrees, which could indicate that the precipitation fell as snow, which is not well captured by the rain gauge at 1.5 m above the soil surface. Some measurements in this period were taken from the nearby experimental station (Geonor, rain gauge, weighing principle) and hence not fully representative of 410 the field. The measured peak in drainage could be caused by snowmelt due to an increase in air temperature. The discrepancy between the estimated and measured could thus be explained by the numerical interpretation of the precipitation by the submodel for snow accumulation and melt incorporated in DAISY. In this sub-model, low temperatures determine the amount of snow-rain mixture and the possible snowmelt is determined by global radiation, air temperature and heat flux on the ground surface (Hansen et al., 1990). The timing of the fast increase in estimated drainage (due to snowmelt) was, therefore, delayed 415 around a week in the simulated results, which may be explained partly by too low a precipitation input to the model.
[ Figure 6 is about here] The other significant disagreement between the model results and the measurements was in mid-December 2006, most likely caused by measurement errors, where measurements indicate zero flux in the drainage, while in this month the model predicted one of the highest fluxes (Fig. 6). In general, however, all six models performed well. The best calibration result was achieved 420 with model with MSET b) + P3, where nMAE was 0.53 for the calibration period for drainage. The comparison of the outcome of applying the two different groundwater tables revealed that all models with P3 outperformed those using P4 as the lower boundary condition in terms of nMAE. It was also clear that P3 set-ups performed much better during the single, sudden flow events such as in October 2006 where the groundwater seems to supply less drainage from the field than with P4 setups. The dynamics of the drainage, including zero fluxes, were perfectly captured applying P3, whereas with P4 the simulated fluxes 425 tended to be high during the period with measured zero flux, probably because the groundwater kept supplying water to the drains. It is also worth pointing out that P4 set-ups performed better for nRMSE and KGE than for nMAE, even though there are visible discrepancies between the observed and modeled flow. This is a clear example of what Muleta (2011) concluded, that using squared deviations tends to overestimate the goodness of fit because peak flows will be given excessive weight. + P3. Additionally, the same identified groundwater-induced drainage surplus appeared here as well with P4, at times when the measurements showed zero flux. Overall, applying the P3 groundwater condition tended to give a better description of the 445 dynamics of the drainage. The reason could be that P3 is installed in a dislocated clay and silt layer that is hydraulically more well-connected to the groundwater aquifer than P4, which is placed in a clay till.

Evaluation of model performance
[ Figure 7 is about here] For situations with fast drainage response to heavy rain, a distinct heavy rain event was screened for model response in SWC at 25 cm depth and drain flow on 8 th June 2003 (Fig. 8). All MSETs and lower boundary setups were able to depict the sudden 450 increase in water content, as well as the immediate water transport to the drain. For the lower boundary, models with P3 were able to simulate the dynamics better. This drainage response is more dependent on the accumulated water above the B-horizon rather than direct transport from the surface. These findings are consistent with the observation of  and Nimmo (2012) that water/solutes can be preferentially transported even though the matrix is only partially saturated below the plow layer during the summer season. Therefore it can be concluded that the DAISY model with Richards' equation,combined 455 with Mollerup's approach (Mollerup, 2010) to PF, is capable of describing PF in unsaturated conditions. Although the model does not account for viscous flow, as a Darcy-based model, it still may indeed represent preferential flow to the drainage over the assessed 10 year period. (Germann, 2018) [ Figure 8 is about here]

Transport pathways 460
The water balance and distribution of different flows and transport pathways to the drain system as a result of the conceptual descriptions of the field were evaluated for the drainage year 2008 -2009. As a percentage of the annual precipitation, the models did not differ significantly from each other (Table 4), although the transport to the drainage with P3 showed higher matrix and macropore water transport to the drain.
[ Table 4 is about here]

465
The optimal fit to the field measurements was achieved by the model with MSET b) + P3 (Fig. 9) estimating a high macropore contribution to drainage (113 mm, 70% of total drainage). In contrast to a) + P3, the matrix contribution to drainage was somewhat higher, although the macropore input remained of the same magnitude (Table 4). Since macropore setting a) does not consider MM macropores, only leaching through the matrix can bypass the tile drain and enter the groundwater. The amount of percolation to the groundwater was similar for all macropores (34-36% of the annual precipitation). This indicates 470 that even though settings b) and c) include the MM and FR macropores as extra transport pathways, no more water percolation to the groundwater (sum of matrix and macropore percolation) appeared. Hence, the simulated degree of contribution from the macropores and matrix to drainage appeared to be a good indication of reality. For setting c) the distribution followed a) and b) when applying P3 and it had a relatively similar outcome to setting b). The implemented fracture gave no extra input to macropore percolation. The simulated results with MSET c) + P4 did not follow the other results.

475
[ Figure 9 is about here]

Bromide transport
Another essential point of the model evaluation was to compare BRtransport of the application 22 nd May 2000 (Fig. 10).
After the application of 30 kg KBr ha −1 ,which is the equivalent of 20.14 kg Br − ha −1 , Fodder beet and Spring Barley were sown in 2000 and 2001, respectively. This showed even though BRis considered a conservative tracer by hydrologists, with 480 limited interactions with the soil microbiome and crop, BRhas been proven to be taken up by crops to a substantial degree (Shtangeeva, 2017). In this modeling case, with P3 as the lower boundary condition, 61%, 57% and 52% of the applied amount was taken up by the sown fodder beet in macropore setting a), b)and c), respectively, whereas the corresponding uptake for P4 was 59%, 60%, and 54%. This amount of BRuptake was not unprecedented in beet crops, as Kohler et al. (2005) found 50% of initially applied BRin harvested sugar beet. However, it must be mentioned that in the present study no crop calibration 485 has been done towards BRuptake.
[ Figure 10 is about here] Two heavy rain events occurred around the day of the application of Br -, the first event causing little drain flow, but the second almost instantly initiating drain flow. The instantaneous occurrence of Brin the drainage at a 5.1 mg l −1 concentration reflects the presence of macropore flow. Figure 10 shows that the Brwas retained throughout the summer period, but started 490 to leach out during the drainage season and the breakthrough ended mid-November. However, due to some retardation in the matrix and in crop residues, later appearances of Brcan be observed. As shown in Figure 10, MSETs with the P3 boundary condition (nMAE BR = 0.85-0.93, nMAE BRC = 0.14-0.21) performed significantly better than any other with P4 (nMAE BR = 1.07-1.19, nMAE BRC = 0.45-0.64), even though P4 set-ups performed marginally better for this specific year in terms of drainage dynamics and cumulative drainage. For model with MSET a) + P4, considerable transport of Brresulted, where 495 K s,B (1.46 cm −h ) was 2.25 times larger than the average K s (0.65 cm −h ) for all set-ups, and therefore more Brleached to the drain through the matrix. For this model, K s,C (0.19 cm −h ) was also one magnitude lower than K s,B , which allowed the water to build up above the C-horizon, which is the reason for the high proportion of matrix water in the drainage (Table 1).
Also, if matrix leaching is included, the different MSETs with P4 resulted in percolation levels that were 1.5-2 times higher than in P3. This supports the hypothesis that the observed lower boundary at P4 had no, or limited, hydrological connection to 500 the field. (Table S2) For all models, the macropore transport of Brfrom the field to the drain substantially outweighed the matrix transport. As suspected, the DM1 macropore class transported the majority of the solute, indicating that the solute did not move further down in the matrix and was being held in the plow layer (A-horizon) to some degree. As Animation S1 shows, after the spraying of the Brthe plume starts to diffuse into the plow layer after the first few hourly rain events, and when the pressure potential 505 exceeds the threshold, all macropores related to the plow layer (DM1,MM1, MM2) are activated and facilitate the transport of Brto the drainage, or below the drain level. On 25 th October 2000, 156 days after the application, the rainy season initiated the Brmovement in the plow layer again and allowed Brto slowly diffuse towards the B-horizon (Animation S2). The animation shows that although the plume had already reached the middle of the B-horizon, the majority of the Brtransport to the drain or below was caused by the remaining concentration in the plow layer. The summarized solute inflow(+) and solute outflow(-) 510 show (Fig. 11) that the transport mechanism was mainly driven by the Brbuild-up above the plow layer, causing large amounts of PF leaching in the soil column.
[ Figure 11 is about here]

Perspective
The present study and results call for further studies to understand the transport of agro-chemicals, such as the optimization of 515 the model with MSET b) + P3 with respect to transport of Br -, nitrate, and applying pesticides. Such future studies could reveal how nitrogen leaching will be affected in the plow layer by mineralization and denitrification processes and further transported to the drainage systems through DM macropores. In the present study, all models proved that the majority of drainage was caused by the DM1 macropores, which collected the water mainly from the plow layer. This reveals that there is a high risk of applied soluble agrochemicals and nutrients, which are incorporated to the plow layer (such as injected slurry), being leached 520 to the tile drain system (Larsson and Jarvis, 1999).

Conclusions
In this study, a one-dimensional model code including different MSETs and a representation of artificial drainpipes in an agricultural clay till field was tested in order to describe the measured water and mass balance, with a focus on drainage. The modeling exercise showed that Darcy flow-based models with coupled dual permeability capabilities could model with high 525 precision for a long-term period while accounting for PF. Two groundwater table measurements, P3 and P4, were collected at the southern and northern borders of the field, respectively, and the results showed that the groundwater table measured at P4 did not truly reflect the water balance of the field, which indicates the importance of a correct lower boundary. The model with MSET b) exposed to groundwater fluctuations measured in the southern part (P3) of the field gave the best description of the drainage of the field. Further, the results of this model study revealed that 70% percent of the overall drainage was supplied 530 via macropores and of the applied tracer, 54% leached directly from the plow layer. Based on the different conceptual settings, the majority of drainage seemed to be primarily the result of rapid precipitation infiltration from the surface to the plow layer, and from there via preferential pulse flow to the drain or below, through macropores. The applied Brtracer test and modeling revealed that in this heterogeneous soil with yearly soil tillage (harrowing and plowing) and, therefore, most likely a temporary increase of the soil porosity in the plow layer, the buried macropores could have a more significant impact on transport than 535 surface-connected macropores.
Video supplement. Animation S1. Simulated bromide transport in the macropore and the matrix media from 22      for all the six model setups (BR = BRdynamics, BRC = Cumulative BRtransport).
37 Figure 11. The summarized Brsolute inflow (+) and solute outflow (-) by depth below the surface, for all six models Tables   825   1 Optimization of hydraulic parameters (vGM) of the matrix for macropore settings (MSETs) a), b) and c). Initial values and selected minimum and maximum uncertainty boundaries are based on laboratory measurements (Table S1). Best calibration for MSET b) shows the result from the best calibration being restricted by the two different fluctuating groundwater tables monitored in wells P3 or P4. The MSET column shows which of the given matrix parameter is active at the given MSET. Sensitivity lists of the MSETs for which a given parameter  Table 1. Optimization of hydraulic parameters (vGM) of the matrix for macropore settings (MSETs) a), b) and c). Initial values and selected minimum and maximum uncertainty boundaries are based on laboratory measurements (Table S1). Best calibration for MSET b) shows the result from the best calibration being restricted by the two different fluctuating groundwater tables monitored in wells P3 or P4. The MSET column shows which of the given matrix parameter is active at the given MSET. Sensitivity lists of the MSETs for which a given parameter is sensitive with given the lower boundaries (P3 and P4). Gray rows are the non-sensitive parameters.  Table 2. Optimization of hydraulic parameters of the macropores for macropore settings (MSETs) a), b) and c). Initial values and selected minimum and maximum uncertainty boundaries are mostly based on Nielsen et al. (2010) (Table S2). Best calibration for MSET b) shows the result from the best calibration being restricted by the two different fluctuating groundwater tables monitored in wells P3 or P4. The MSET column shows which of the given macropore parameters is active at the given MSET. Sensitivity lists of the MSETs for which a given parameter is sensitive with the lower boundaries (P3 and P4) given. Gray rows are the non-sensitive parameters.   Table 4. Accumulated water contribution to drainage from matrix and macropores (DM1, DM2 and the sum of the two) and percolation through the matrix and macropore at 1.1 m depth (depth of drain) in mm, and as a percentage of the annual precipitation (997 mm) for the period 2008-2009 for all six models (MSET a), b) and c) applying groundwater tables from P3 and P4).