Three Two-Dimensional Approaches for Simulating the Water Flow Dynamics in a Heterogeneous Tile-Drained Agricultural Field in Denmark

across hydro-topographical subcatchment selected study area spatial the geological, stratified heterogeneity. Three two-dimensional models incorporating preferential water flow used to examine the ability of each approach to predict the drainage dynamics and water balance in a clayey subsurface drained agricultural field: (i) a single-porosity model, (ii) a dual-porosity model, and (iii) a dual-permeability model. The initial parameterization of the models relied on in situ obser-vations of soil structure, soil hydraulic property maps of Denmark and pedotransfer functions. The models, once calibrated against observed drainage data, were used for the reproduction of the subsurface drainage water dynamics for different time periods for the specific agricultural landscape. The different conceptual approaches gave a fairly good fit to the hydrograph features for both the calibration and validation datasets. The predictive ability of all models for the calibration dataset was fairly similar, while for the first and second validation datasets better results were obtained with the dual-porosity and dual-permeability models, respectively. The calibrated values of the effective parameters had been examined using laboratory data and in situ measurements of the groundwater level to evaluate how well these parameters physically represent the flow domain.

N on-equilibrium preferential flow models are increasingly used in modeling the water flow dynamics and the nitrate fate across tile-drained agricultural fields (Mohanty et al., 1998;Larsson and Jarvis, 1999a;Šimůnek et al., 2003;Gärdenäs et al., 2006). Tile drainage is essential for the efficiency of agricultural production, since it removes excess soil water and protects the crops from anoxic conditions (Fausey et al., 1987). However, tile drains have a significant impact on the hydrology and flow patterns in a drained agricultural field, operating as conduits and constituting the link between the soil profile and surface and/ or subsurface water receptors (Frey et al., 2016). Several studies have examined how surface-applied agricultural fertilizers can be transported through the vadose zone, increasing the risk of surface water bodies, and groundwater contamination (Larsson and Jarvis, 1999b;Jørgensen et al., 2004;Gärdenäs et al., 2006). Preferential flow is an important process in the transport of water and contaminants in the soil. The flow process in clay or structured loam soils occurs to a large extent in macropores such as biopores (e.g., earthworm burrows and plant roots) or in planar fissures formed by shrinkage and swelling (Jarvis et al., 2009;Iversen et al., 2012;Larsbo et al., 2014). In addition, the flow process occurs in interaggregate packing voids created by topsoil management (Beven and Germann, 1982;Jarvis, 2007). Gerke and van Genuchten (1993a) found that these microscopic structures influence the water movement and solute transport by developing non-equilibrium flow fields with different magnitudes of velocities. The definition of a macropore is equivocal (Beven and Germann, 1982) and a wide range of equivalent cylindrical diameters have been proposed as a definition for macropores. Luxmoore et al. (1990) stated that: "The term macropore includes all pores in a profile that are (generally) drained at field capacity, with the latter being 1 mm or more in equivalent diameter". Jarvis (2007) defined macropores as pores of equivalent cylindrical diameters larger than ~0.3 to 0.5 mm which drain at a suction level of -10 cm H 2 O. In addition, Moldrup et al. (2003) defined macropores as pores with an equivalent pore diameter of more than 30 mm which drain at a suction level of -100 cm H 2 O. Moreover, several studies have defined macropores as pores larger than 50 mm (Carter, 1988;Pagliai et al., 2004;Passoni et al., 2015). During the last decades, a variety of sophisticated models have been developed for describing the macropore and non-equilibrium flow at different scales, for exapmle, dual-porosity, dualpermeability, multi-porosity, and/or multi-permeability models (Šimůnek and van Genuchten, 2008).
Nitrogen transport through the subsurface of the soil is a complex process. While infiltrating into the soil, ammonium and nitrate can be taken up by the plants and ammonium can be transformed to nitrate by nitrification processes. Jørgensen et al. (2004) found transport of nitrate to be controlled by flow via macropores in clay soils, and Mohanty et al. (1998) and Larsson and Jarvis (1999a) showed the importance of macropore flow in relation to nitrate transport in structured soils. They discovered that the presence of macropores and their direct connection to the tile drainage network are strongly linked to nitrate contamination problems. In contrast, Bergström (1995) concluded from leaching experiments that nitrate is better distributed in the soil matrix and water rapidly moving through macropores will not carry significant amounts of nitrate compared to pesticides. Lee et al. (2014) investigated the effect of soil texture on the transport of nitrate N. In that study, the breakthrough curves of nitrate N for a coarse-texture soil followed a typical Gaussian distribution, indicating matrix-dominated transport, while the breakthrough curve for a structured soil revealed the contribution of preferential flow to solute transport. Therefore, since redox transformation processes are closely linked to the spatially distributed water content and the residence time of nitrate in the soil, the process of nitrogen reduction and evolution of contaminant transport can be simultaneously controlled by the flow via the soil matrix and the preferential flow pathways because of subsurface structural heterogeneity.
The total land area in Denmark under agricultural use is approximately 60% with more than 90% as arable farming (Kronvang et al., 2017). It is assumed that approximately 50% of the arable land is drained (Olesen, 2009) and therefore associated with nutrient contamination of freshwater systems via drains (Billy et al., 2011). Simultaneously, part of the percolated nitrate N enriched water bypasses the subsurface (tile) drains and recharges local aquifers by following natural flow pathways and finally, discharges in surface water systems (Hinsby et al., 2008). Consequently, increased nitrate concentrations are measured in Danish groundwater, streams, lakes, and fjords. Most of this nitrate is believed to originate from nonpoint pollution from agricultural areas due to the major loss of N from cropping systems (Askegaard et al., 2011). A major challenge for Danish water resource management is therefore to reduce the nitrate load to lakes and coastal waters from agricultural areas. Agricultural practices have been strictly regulated in Denmark since the 1980s to reduce nitrate leaching (Hansen et al., 2014). In 2013, more targeted regulations were deemed necessary, in which the local attenuation capacity for nutrients in the landscape would be taken into account (Kronvang et al., 2017). Therefore, a locally differentiated future regulation taking into account the spatial variability of the natural transformation of nitrate in the subsurface has been set as future goal (Hansen et al., 2014). This requires better water resource management tools, regarding the modifications of the currently performed agricultural practices.
Modeling seasonal water flow and nitrogen reduction across hydro-topographical gradients and tile-drainage sequences in agricultural landscapes requires consideration of all the complex hydrological and biogeochemical dynamic processes that take place in the vadose as well as the saturated zone. However, before considering chemical species, an accurate understanding of the water balance is essential.
Depending on the soil processes that need to be modeled, selection of the most suitable approach for generating a hydrological model is essential. The objective of this study was to examine three different model approaches in terms of their capabilities and weaknesses for simulating the soil water regime and water balance across hydro-topographical gradients and drainage sequences in a Danish catchment. Three models in HYDRUS-2D were used for that purpose which all incorporated preferential water flow.

Study Area
The selected drainage catchment (10.4 ha) is located within the Norsminde catchment (100 km 2 ) in eastern Jutland, Denmark (55°59¢ N, 10°4¢ E; Fig. 1). It is characterized by a rolling hilly terrain where the surface elevation varies from 95 to 72 m above mean sea level. Clayey till facies are located down to 3-m depth below which are clayey and sandy glacial deposits De Schepper et al., 2017). The deeper geology is characterized by Miocene deposits on top of Paleogene fine-grained marine marls and clays (De Schepper et al., 2017).
An electromagnetic survey was conducted in the study area using a DUALEM-21S sensor (DualEM, Milton, ON, Canada). Based on the soil maps produced (Fig. 2), the optimum locations considering the spatial heterogeneity of soil properties were se- lected for installation of 54 partially penetrating piezometers at different depths to monitor the piezometric head of the more conductive layers. Additionally, the field was equipped with three fully penetrating piezometers (approximately 2-m depth, hereafter referred to as P17, P18, and P19) for monitoring the groundwater level once per month. During the installation of the piezometers, 76 soil samples were collected from various depths, 22 of which from two depression zones. The depression zones were identified at the studied field from the Digital Elevation Model (DEM), DUALEM measurements and in situ observations. Sixteen representative soil samples were used for soil texture and soil organic matter analysis. For the determination of soil texture, a combination of sieving and the hydrometer method was used (Gee and Or, 2002). Total organic C was determined on a Leco Carbon Analyzer coupled to an infrared CO 2 Flash 2000 NC detector (Thermo Fisher Scientific Inc., MA). Based on the soil texture analysis, the drainage subcatchment was classified as mainly sandy loam, except for a few subsoil samples that ranged from sandy loam to clayey loam (Table 1; Fig. 3).
This difference in soil texture with depth reflects the pedological clay illuviation processes characterizing this clay moraine catchment. In addition, the hydro-topographical gradients resulted in a redistribution of soil particles and soil organic matter and thus has led to higher clay and organic matter contents in the two depression zones. The borehole descriptions during the installation of the piezometers were used for defining the average thickness of each soil horizon (Fig. 4), which was also used to set up the conceptual model for the studied field. Furthermore, a total of 36 undisturbed 100-cm 3 (3.5-cm height and 6.1-cm diameter) soil cores were sampled from 0-to 20-cm depth at different locations in the drainage subcatchment ( Fig. 2) for soil-water retention analysis. For the soil water retention (SWR) curves from -10 to -1000 cm H 2 O, the soil cores were slowly saturated with tap water from the bottom and subsequently drained to the intended matric potentials utilizing a sandbox (-30 and -100 cm H 2 O), hanging water column (-300 cm H 2 O) and the Richards pressure plate apparatus (-1000 cm H 2 O). At each measurement step, the volumetric water content (q, cm 3 cm -3 ) was calculated by multiplying the gravimetric water content with soil bulk density and dividing by the density of the tap water. Bulk density (g cm -3 ) was determined for the 100-cm 3 core samples by oven-drying at 105°C for 24 h.
The drainage catchment is tile-drained (10-cm diameter) at a mean depth of 1.2 m. The cone of depression is assumed to be 10 m from the center of tiles equal to half the distance between the mapped tile drains. Drainage discharge is monitored by an    Event101A data logger (MadgeTech, Warner, NH) recording the pulse every 100 L by a bidirectional, battery-operated, electromagnetic Waterflux 3070 flow meter (KHRONE, Duisburg, Germany) connected to the drainage network outlet and is continuously converted into hourly data.

Flow Models Description
The modified HYDRUS-2D finite element code (Šimůnek and van Genuchten, 2008;Šimůnek et al., 2016) was used for the simulations of the tile drainage dynamics. In this study, one equilibrium and two non-equilibrium modeling approaches were used: a single porosity model using modified soil hydraulic properties to account also for the contribution by macropores (Vogel and Cislerova, 1988), a dual-porosity model (Šimůnek and van Genuchten, 2008;Šimůnek et al., 2016) and a dual-permeability model van Genuchten, 1993a, 1996). All these approaches describe the variably saturated water flow through the porous system using the Richards' equation (Šimůnek and van Genuchten, 2008). The main differences between the models are in their conceptualization of the pore system and the part of the porous medium in which water flow occurs. The approaches were selected to be two-dimensional since, the two-dimensional approaches are suitable for studying tile-drained fields, in which the tile-drains are ascribed by their diameter, depth and spacing (De Schepper et al., 2017).

Single Porosity Approach Incorporating the Contribution of Macropore Flow
In the single porosity model, the variably saturated flow is described using Richards' equation: where q (L 3 L -3 ) is the volumetric water content, t (T) is the time, h (L) is the suction, S (T -1 ) is a sink term accounting for root water uptake, x i (i = 1,2; L) are the spatial coordinates and K ij are components of a dimensionless anisotropy tensor K A . The K (L T -1 ) is the unsaturated hydraulic conductivity function. The soil water retention, q(h), and the hydraulic conductivity, K(h) for variably saturated conditions are characterized using the van Genuchten-Mualem model (van Genuchten, 1980): where q r and q s (L 3 L -3 ) are the residual and saturated water contents, a (L -1 ) is related to the inverse of the air entry suction, h (L) is the suction, n (-) is a measure of the pore-size distribution, m = 1-1/n (-), K s (L T -1 ) is the saturated hydraulic conductivity, l (-) is the tortuosity parameter, and S e (-) is the effective saturation degree given by: The contribution of the preferential flow patterns to the total water flow has been taken into account using the model suggested by Vogel and Cislerova (1988). The hydraulic conductivity function is here modified to add flexibility in the description of the hydraulic properties near saturation (Šimůnek and van Genuchten, 2008). The modified van Genuchten model assumes that unsaturated hydraulic conductivity for a capillary dominated flow can be described by Eq.
[5] until the suction reaches a critical value h k and the hydraulic conductivity is equal to K k = K(h k ). Afterward, a linear increase of hydraulic conductivity from K k to K s for suction h k < h < h s is assumed to account for non-capillary, macropore-dominated flow (representing effective macroporosity or other structural features) until it reaches a suction equal to h s and hydraulic conductivity K s .
This approach is also in accordance with Jarvis (2007), Børgesen et al. (2006) and Iversen et al. (2011), who studied the link between the hydraulic conductivity at saturated and near-saturated status of the soil and the potential of water flow through soil macropores.

Dual-Porosity Approach
For the dual-porosity approach, the porous medium consists of two interacting domains. One linked with the macroporous system and one to the micropores or the intra-aggregate pores (Šimůnek and van Genuchten, 2008;Šimůnek et al., 2016). Therefore, the dual-porosity model conceptualizes the flow domain into a mobile (subscript "f ") and an immobile (subscript "m") region (Šimůnek and van Genuchten, 2008). In HYDRUS-2D, the domains are not geometrically separated, but are considered as two overlapped regions in the same soil volume (Köhne et al., 2009). The mobile region (the macropores; q f [L 3 L -3 ]), is assumed to be the part of the soil where water is restricted to flow, while the immobile domain (the soil matrix; q m [L 3 L -3 ]) has zero saturated hydraulic conductivity (stagnant), allowing only drying or rewetting processes. Thus, the total water content is given by: Water flow in macropores is described using Richards' equation and a mass balance equation to describe moisture dynamics in the matrix (Šimůnek and van Genuchten, 2008): where S f and S m (T -1 ) are sink terms for both regions, and G w is the exchange rate (T -1 ) for water from the macropores to the matrix and vice versa. Equation [9] denotes that the water transfer term is proportional to the difference in effective saturation, S e (0 ≤ S e ≤ 1), between the macropores and the matrix (Kohne et al., 2004) and w (T -1 ), which is a first-order mass transfer coefficient: where q f and q m (L 3 L -3 ) are the mobile and immobile water contents with residual and saturated constraints of q r,f , q r,m , q s,f , and q s,m , respectively. For the calculation of the effective saturations, it is assumed that the retention curves of the mobile and immobile domains are similar (Kohne et al., 2004).

Dual-Permeability Approach
Similar to the dual-porosity model, the dual-permeability model considers the total pore space as a two-region system (macropores or inter-aggregate pores and micropores or intraaggregate pores). In the dual-permeability model there is no limitation having direct flow of the liquid phase between the aggregates themselves in contrast to the dual-porosity model where this is not possible (Šimůnek and van Genuchten, 2008).
Richards' equation is applied to both domains. The flow equations for the macropore (subscript "f ") and matrix (subscript "m") pore systems are given by van Genuchten, 1993a, 1993b): where w ( -) is the ratio of the volumes of the macropore or fracture domain and the total soil system. It is important to mention that water contents q f and q m in Eq.
[11] and [12] refer to the local water contents of the two regions, while in the dual-porosity model (q F , q M in Eq. [13]) they refer to the total pore space (Šimůnek et al., 2012).
Equation [14] indicates that the water transfer term in the dual-permeability approach is proportional to the suction difference between the macropores and the matrix (Gerke and van Genuchten, 1993a) and the a w (L -1 T -1 ), which is a first-order mass transfer coefficient. For this approach, the water release characteristics of the two domains are represented by different soil water retention curves for each of the defined regions (Kohne et al., 2004): where b (-) is a shape factor that depends on the geometry, d (L) is an effective diffusion path length (i.e., half the aggregate width or half the fracture spacing), and g w (-) is a scaling factor. In addition, K Int (L T -1 ) is the hydraulic conductivity of the macropore-matrix interface, being equal to the arithmetic average of the fracture-matrix interface hydraulic conductivities at a given h f and h m in the two domains (Gerke and van Genuchten, 1993b):

Flow Domain, Boundary Conditions, and Initial Conditions
For all three modeling approaches, the simulated domain had a width of 20 m and a depth of 1.3 m. Furthermore, one drain at 1.2-m depth was placed in the middle of the flow domain (Fig.  4). The selection of the specific conceptual layered domain (Fig.  4) relied on the assumption of a 10-m cone of depression and the presence of a very low-permeability layer below 1.3-m depth. The latest was based on the DUALEM measurements, the soil sampling during the installation of the piezometers as well as the borehole descriptions. In addition to the previous, the soil texture analysis was also used for dividing the soil profile into four distinct horizons (Fig. 4). The top two horizons (Ap and E) had thicknesses of 0.25 m, followed by two clay accumulation horizons with pseudogley characteristics Btg1 and Btg2, both having thicknesses of 0.4 m. The relatively impermeable layer under a calcareous boundary was designated a Cr horizon.
Boundary conditions for the simulations of the three approaches included no-flux boundary conditions at the bottom and lateral boundaries of the system. However, the assumption of a no-flow Neumann boundary condition at the bottom is expected to influence the accuracy of the simulated cumulative drainage discharge, since a very small part of the infiltrated water would have bypassed it. Atmospheric conditions at soil surface were used (hourly precipitation, evaporation and transpiration rates). The tile drain was treated as a seepage face boundary condition considering that the suction along the drain was equal to zero when the surrounded soil was saturated. In contrast, when the water content of the surrounded soil was at a suction below zero, the tile drain was treated as a no flow boundary condition (Šimůnek and van Genuchten, 2008;Šimůnek et al., 2016). The simulated domain was represented by 10,491 triangular finite elements. Close to the surface and around tile drains the nodes were given higher nodal densities to achieve a mass balance error of less than 1% for water flow during the precipitation events.
The same initial conditions were set for all models and the determination of the water level at t = 0 relied on the measurements from piezometers P17, P18, and P19. The DUALEM data and the DEM were used to calculate three weighted factors for estimating the initial average water level for the tile drain catchment using the measurements of groundwater level from P17, P18, and P19. This was based on the soil texture, the topographical gradients, and the location of each of the piezometers. Thus, the groundwater level at the beginning of the simulation was set to be 10 cm above the bottom boundary and hydrostatic equilibrium was applied within the domain.

Drainage Discharge, Meteorological Conditions, and Root Water Uptake
The available data was subdivided into calibration and validation datasets to determine and test the effective parameters for each of the three selected model approaches. For both selected time periods, there was not surface runoff based on visual observations.

Drainage Discharge Calibration-Validation Datasets
The selected period for the calibration of the models was from 3 May 2013 to 18 June 2014 (9744 hourly values) as the most continuous data series of tile-drainage discharge for the tile drain catchment including wet and dry periods.
For the validation of the models, two periods were selected based on the continuity of measured tile drainage discharge data and the desire to include for both simulation periods a wetting and a drying period. Especially the first selected dataset for the validation period is characterized as a period with intensive drainage flow and with a wide range of fluctuations in the discharge. The second dataset was selected to test the simulation capability of the model of the transition from very dry to wet conditions. The first validation period was from 13 Nov. 2014 to 13 May 2015 (4367 hourly values) and the second validation period was from 4 Aug. 2016 to 25 Oct. 2016 (1992 hourly values).

Input Variables
Hourly precipitation data was used from two nearby meteorological stations located respectively 3 and 6 km from the studied area. Potential evapotranspiration was calculated with the Penman-Monteith equation, using as input variables data extracted from the meteorological station in Foulum, Denmark (56°29¢ N, 9°35¢ E) (distance 62.5 km). The daily potential evaporation and potential transpiration were estimated from the water balance model EVACROP (Olesen and Heidmann, 1990). This is necessary as the HYDRUS-2D has no crop-growth module and cannot subdivide the potential evapotranspiration data into actual water balance components. For the conversion of the daily potential evapotranspiration into hourly data, a Gaussian distribution was used and the calculated potential evapotranspiration data were normally distributed between the hours of 5:00 AM and 8:00 PM (Bakhtiari et al., 2006) with maximum density between 12:00 PM and 1:00 PM. Winter wheat (Triticum aestivum L.) was used as a model crop for the whole simulation period and the sowing date was always considered to be the first day of September due to unavailability of actual management data.
Root water uptake was simulated using the reduction model of Feddes et al. (1978) and the critical stress index was set to 0.75 to take into account the compensated root water uptake. The maximum rooting depth and the depth of maximum root density were set to 1.2 m (Palosuo et al., 2011) and 0.7 m (Gärdenäs et al., 2006), respectively. The threshold values of the Feddes et al. (1978) reduction function were based on the default values for wheat in HYDRUS-2D (Wesseling, 1991) and afterward, they were adjusted for winter wheat according to Gärdenäs et al. (2006). Therefore, the threshold values for anoxic conditions (h 1 and h opt ) were set to 0 and -0.1 m, h 3 to -15 m, while the value for the parameter h 4 was set to -160 m. For both calibration and validation of the models, the hourly potential evaporation and potential transpiration were estimated similarly.

Effective Properties Parameterization
The initial parameterization for the effective properties relied on field measurements at the drainage catchment, on previous studies, and on a literature review (Table 2). For this study, the macropores are defined according to Carter (1988); Pagliai et al. (2004); Passoni et al. (2015). Initially the dual porosity model was calibrated using as an objective function the observed drainage discharge. Then, the effective parameters from the dual-porosity model were inserted into the dual-permeability model as fixed values and the remaining parameters for the matrix domain were calibrated. Finally, the relevant parameters from the calibrated dual-permeability model were inserted into the single-porosity model, which was calibrated by varying K s . This workflow is illustrated in Fig. 5.

Dual-Porosity Model
For the dual-porosity model, the total porosity (f) for the Ap horizon was estimated from the measured dry bulk density (r b = 1.49 g cm -3 ) and the equation f = 1r b /r d (Brady and Weil, 2013), where r d is the particle density and was assumed to be equal to 2.65 g cm -3 (McBride et al., 2012;Schjønning et al., 2017). For the other three layers, the r b was set according to Børgesen and Schaap (2005). Saturated water content, q s , was then assumed to be equal to the total porosity of the porous medium. Total porosity was additionally evaluated whether it presented a decreasing trend through the soil profile, similar to other studies in sandy loam fields in Denmark (Christiansen et al., 2004). The q f was estimated from a pedotransfer function (PTF) (Stolf et al., 2011) and from a literature review (Katuwal et al., 2018) for evaluation of the macroporosity. Stolf et al. (2011) developed a PTF for estimating macroporosity (defined as the Table 2. Initial values of effective parameters for all the geological layers in 2D modeling approaches: Single porosity (SP), dualporosity (DP), and dual-permeability (DUP) for four depths.   pores with an equivalent pore diameter >50 mm) as a function of sand content, soil bulk density, and particle soil density. Katuwal et al. (2018), on the other hand, measured a macroporosity of 0.086 cm 3 cm -3 using x-ray computed tomography in a topsoil for a similar field as the studied field in terms of soil texture and crop management (Estrup, Denmark) (Paradelo et al., 2015) where macropores were defined as having an equivalent pore diameter of more than 1.2 mm. A significant decrease of macroporosity with depth takes place in the plowed Ap horizon (0-25 cm) or through the combined Ap and E horizons (Asare et al., 1999). For this study, the largest difference between the macroporosity across the soil profile was set from the A p to the E horizon. For the van Genuchten parameters a and n, the values from Gärdenäs et al. (2006) were used as initial parameterization, and the parameter l was set to -1 for all the approaches according to Leij et al. (1996). Iversen et al. (2011) estimated the potential of water flow through soil macropores. These authors developed PTFs for estimating the hydraulic conductivity of macropores, K s,f , using the difference of the predicted soil hydraulic conductivity at a matric potential of h = 0 and h = -10 cm H 2 O [K s and k(-10 cm H 2 O)] for a variety of different Danish soil types.
In the present study, it is assumed that the decrease in hydraulic conductivities from saturation to h = -60 cm H 2 O can be regarded as equal to the decrease in hydraulic conductivity from saturation to h = -10 cm H 2 O considering that at those suction levels the macropores had been drained and do not contribute to the water transport. Therefore, values of k(-10 cm H 2 O) based on the PTF of Iversen et al. (2011) were used as an estimate of k(-60 cm H 2 O) and their difference from K s used as initial K sf input for each layer. Moreover, the soil hydraulic properties were assumed homogeneous and isotropic in each layer and, therefore, K xx and K zz were both equal to the unsaturated hydraulic conductivity function, K(h).

Dual-Permeability Model
For the dual-permeability approach, the hydraulic properties for the soil matrix were estimated, while for the macropores the calibrated effective parameters from the dual-porosity model were applied. For the ratio w of the volumes of the macropore domain and the total soil system, the incipient value was selected according to Gärdenäs et al. (2006) in order for the derived value for the A p macropore-saturated water content to be equal to 0.75. The two van Genuchten parameters (a, n) for the soil matrix were initially defined according to Børgesen and Schaap (2005). In that study, parametric PTFs for the estimation of van Genuchten (1980) water retention parameters were developed using the Danish profile database of soil hydraulic properties and utilizing as predictors the soil texture, bulk density, soil organic matter content, and horizon information.
The initial estimation of the saturated hydraulic conductivity of the soil matrix (K s,m ) for the four defined layers in the dual-permeability approach was set according to Nielsen et al. (2018), in which saturated hydraulic conductivity prediction models for top-and subsoil were developed. In that study, PTFs were derived using K s measurements of 100-cm 3 intact soil cores from arable land across Denmark. The initial estimated value of K s,m for the topsoil was fairly compared with measurements from Estrup (Karup et al., 2016) using only soil cores showing signs of a low biopore activity to minimize the effect of the macropores. According to Messing and Jarvis (1995) the undisturbed small core samples (10-cm height and 7.2-cm diameter) underestimated the hydraulic conductivity in the suction range -15 cm H 2 O to saturation compared to larger samples (approximately 45-cm height and 30-cm diameter) because macropores appeared less frequently and might have been underrepresented in these samples. Furthermore, the hydraulic conductivity decreases by one to three magnitudes for a small suction increase from h = 0 to -10 cm H 2 O, owing to the effects of structural macropores (Clothier and Smettem, 1990;Messing and Jarvis, 1995). Iversen et al. (2001) also found a relatively large scaling exponent between small (100 cm 3 ) and large (6280 cm 3 ) soil samples for the structured soils at Estrup. In the study by Nielsen et al. (2018), even smaller soil cores (3.5-cm height and 6.1-cm diameter) from agricultural landscapes across Denmark were used, and therefore it was assumed that a satisfactory initial estimation for the K s,m at -60 cm H 2 O can be ascribed. The a w , b, and g parameters were set to 3 and 0.4, respectively (Gerke and van Genuchten, 1993a;Gärdenäs et al., 2006;Šimůnek et al., 2012). The saturated hydraulic conductivity of the macropore-matrix interface, K s,Int , for all horizons was assumed to be 1% of the saturated conductivity of the matrix van Genuchten, 1993a, 1993b).

Single-Porosity Model
The single-porosity model included the use of modified soil hydraulic properties to account for increased flow near saturation. The initial parameterization relied on the results from the calibration related to the matrix hydraulic properties of the dual-permeability model. The derived values from the calibration of the dual-permeability model were validated with the measured soil water retention data, assuming that the hydraulic properties of the extracted core samples were characterized by a single flow domain.
The threshold water content (q k ) determining the contribution or the non-activation of the macropores was defined based on the soil water retention data (derived from the calibration of the dual-permeability model) at a suction, h k , of -60 cm H 2 O. This was combined with the equation e 60 = fq -60 cm H2O (Resurreccion et al., 2008), where e 60 corresponds to the air content at -60 cm H 2 O of suction and f to the total soil porosity. The hydraulic conductivity K k for h = -60 cm H 2 O was set equal to the K s,m as it was finally calibrated in the dual-permeability model, while the incipient value for K s was estimated according to Iversen et al. (2011).
The final parameterization was achieved by repeated stepwise manual calibration utilizing the observed drainage discharge as the objective function and sustaining the input parameters within a physical range that would also correspond to the analyses such as soil texture and SWR analysis. An additional objective function was the non-appearance of surface runoff during the entire simulation. The outcome from each calibration was evaluated using the R 2 and the RMSE. Prior to that, a sensitivity analysis was performed to identify the important parameters for each of the approaches.

Statistical Evaluation of the Models
For the statistical evaluation of the models, the coefficient of determination (R 2 ), root mean square error (RMSE), the ratio of performance to inter-quartile distance (RPIQ) and the Akaike's information criterion (AIC) were used. The RMSE is a measure of the differences between predicted values and observed data, defined as: where N is the number of samples,  i x is the predicted values and x i is the observed data.
The RPIQ represents the spread of the population regardless of the distribution (Bellon-Maurel et al., 2010). The index is based on quartiles where Q 1 is the median of the lower half of the data set and Q 3 is the median of the upper half of the data set: where IQ is the inter-quartile distance and gives the range of dispersion around the median, and SEP is the standard error of prediction. The RMSE is used instead of the SEP.
To account for the number of model parameters when comparing model performance, Akaike's information criterion (AIC) was used (Akaike, 1973;Carrera and Neuman, 1986;Hwang et al., 2002).
where d i is the difference between simulated and observed drainage discharge, and k is the number of model parameters. Specifically, only the number of calibrated model parameters in each approach had been used considering the overlap parameters among the models. Smaller AIC values indicate better model performance (Minasny et al., 1999).

Measured and Simulated Drainage Discharge Dynamics
The measured drainage discharge was characterized by a large number of distinct peaks and short recessions indicating rapid water flow through macropores. Table 3 presents the cali-brated effective parameters. All of the used approaches were able to simulate the main attributes of the flow hydrographs with the desired degree of accuracy and the model predictions broadly matched the timing and general shape of the observed water outflow. However, the output from the different approaches showed similar difficulties in simulating the same specific periods of time (Fig. 6a-c) where all the models captured the distinct, high peak flows but simulations consistently failed to match the low continuous discharge rates during the dry seasons.
Initially, from 15 to 17 May 2013, the three simulations predicted similar lower values compared to the observed discharge rates, which was approximately 2 L s -1 . The highest degree of underestimation was observed for the single-porosity model. This difference between the observed and modeled drainage flow can be caused by a discrepancy between the local rainfall and the input values from the nearby meteorological stations. Another explanation is uncertainty in the initial conditions. All three approaches matched well the peak (~34 L s -1 ), from 22 to 30 May 2013. Especially the dual-porosity model showed fairly good agreement with the observed drainage, with a prediction error of 0.6%. Subsequently, the models failed to predict a steady base outflow with relatively small fluctuations (~0.2 L s -1 ) during the dry period of the simulation. This failure may relate to an unknown source that recharges continuously at a point in the drainage network. Such unknown source could be explained by high groundwater levels present in the two defined depression zones, which may contribute to the continuous, low drainage discharge. The period from 29 Oct. 2013 to 8 Mar. 2014 was characterized by intensive drainage output showing wide range of fluctuations in discharge. The dual-permeability model showed here the best predictive ability of the drainage dynamics, while the dual and single-porosity models slightly overestimated the drainage outflow. From 29 Oct. 2013 to 5 Dec. 2013, the drainage discharge was poorly predicted by all models. The reason may be an overestimation of precipitation from the nearby meteorological stations, and maybe also an underestimation of evapotranspiration due to the limited knowledge of crop management in the studied area.
For all models an overestimation of the cumulative drainage output was observed (Fig. 6d), with similar cumulative bias (Table 4) compared to the observed drainage discharge. This could be explained by the assumption of no flux at the bottom boundary of the simulated flow domain. Furthermore, the drainage catchment was in the modeling defined as a systematically drained agricultural field by using a simplified conceptual domain for modeling the entire drainage subcatchment. However, the field is only partly tile-drained (Fig. 1), which may have caused recharging of the local subjacent aquifer bypassing the tile-drainage network De Schepper et al., 2017).
Generally, the predictive ability of the three approaches in relation to the timing and shape of the outflow hydrographs was good. This is confirmed by the range of the coefficient of determination (R 2 ) being similar for the three models. The performance statistics of all models for the calibration period Table 3. Calibrated effective parameters for the matrix (subscript "m") and macropore (subscript "f") regions of the different soil materials as used in the various modeling approaches: Single-porosity, dual-porosity, and dual-permeability models.

Single porosity
Depth q r q s a n K s q r,f q s,f a n K s,f q r,m q s,m w are presented in Table 4. For the three models, R 2 was equal to 0.84, while the RMSE for the single-porosity approach was slightly higher (~5%) than for the dual-porosity and dualpermeability models. Based on the RPIQ index, a statistically better performance was seen for the dual-porosity model with value of 2.27 compared with 2.24, and 2.12, respectively, for the single-porosity and dual-permeability models. However, the superiority of the dual-permeability model is evident from its ability to match peak discharge values and recession limbs more precisely and its smoother transition from wet to dry conditions. This is probably because of its inclusion of interaggregate water-exchange, which was not included in the dualporosity model, resulting in an abrupt transition. Moreover, the different SWR curves for the mobile and immobile regions allowed the stepwise behavior of the dual-permeability approach, while for the dual-porosity the similarity in the SWR curves of these two different domains prevented it. Finally, the AIC values indicated that the dual-permeability model was the best performing model, where for the single-porosity approach the k value was 4 and, respectively, 28 and 17 for the dualporosity and dual-permeability approaches.

Model Validation Simulation of Discharge Dynamics for Different Time Periods
The validation of the three calibrated approaches using data from two different time periods was in agreement with the observed measurements. Table 4 shows the performance statistics for the approaches used. The tile drainage output for the first validation dataset (Fig. 7) was simulated by the dual-porosity model with an R 2 value of 0.77 and RMSE of 1.62 L s -1 . For the dual-permeability and single-porosity models these values were 0.71, 0.75 and 1.67, 1.96 L s -1 , respectively. The AIC values were 1.67 × 10 4 , 1.69 × 10 4 and 1.83 × 10 4 , respectively, for the dualporosity, single-porosity, and dual-permeability models.
For the second validation period, the models also showed fairly good predictive ability (Fig. 8). However, higher discrepancies were observed between the three approaches. For the dual-porosity model drainage discharge simulation achieved an R 2 value of 0.83 and RMSE of 1.59 L s -1 , while for the single and dual-permeability models these values were to 0.77, 0.85 and 2.02, 1.46 L s -1 , respectively. The AIC values for the dual-porosity, single-porosity, and dual-permeability models were 7.56 × 10 3 , 8.45 × 10 3 and 7.19 × 10 3 , respectively.
The RPIQ index for the first validation period indicated that the single-porosity and dual-porosity models more accurately predicted the drainage discharge, while for the second validation the dual-permeability model showed a better performance. These results were expected, since the first validation was a time period with intensive tile drainage discharge and the dual-porosity and single-porosity approaches were capable of producing high peaks and simulating water flow in the fastflowing region. In contrast, the dual-permeability model underestimated these high peaks, probably due to the inclusion of water flow in both regions, thereby reducing peak flow rates. Nevertheless, the dual-permeability model had a higher RPIQ value for the second validation compared to the other two models, indicating that it more accurately simulated the transition from dry to wet conditions.

Validation of Soil Hydraulic Properties with Field and Laboratory Measurements
A successful fitting of the hydrographs for both calibration and validation datasets does not necessarily ensure the representation of the actual internal flow and transport pathways within the studied drainage subcatchment. Hence, the calibrated values of the effective parameters (Table 3) required validation to evaluate how well these parameters physically represented the flow domain. First, we evaluated the van Genuchten parameters (a, n) for the macropores (>50 mm) and the matrix (<50 mm) for each of the four layers in the model. For the macropore characteristics, the calibrated values are comparable to coarse sand. The selected values for the van Genuchten parameters conferred a flatter shape to the soil-water retention curve for the macropores (Fig. 9a). For parameter a, the selected values allowed the macropore domain to be drained or become fully saturated at very low suction. In addition, the high values for parameter n trigger a sharp decrease in water content for very low suctions in the fast-flow region (Cheng et al., 2017). This behavior compares well with the given definition of the macropores for this study, since the macropores were almost totally drained at -60 cm H 2 O. Thus, the hydraulic conductivity of the macropores was allowed to decrease from an extremely large saturated value to nearly zero at small tensions. Simultaneously, the slight decrease of parameters a and n with depth reflect the clay illuviation in these soil profiles. Figure 9b shows the SWR curve for the matrix domain, and its shape agrees with the assumption that pores with an equivalent pore diameter of less than 50 mm start becoming drained at a suction level of -60 cm H 2 O.
The water release characteristics of the soil as a single domain were validated using the soil water retention data. The calibrated van Genuchten parameters (a, n) and the predicted q s and q r were inserted into the van Genuchten function to predict the SWR curve for each layer in the model and then compared with the measured q(h). Taking into account that soil samples for soil water retention measurements were collected from the topsoil only, the van Genuchten parameters were statistically examined only for the first geological layer in the models (0-25 cm). For the other three layers, it was evaluated only if they imitate the hydraulic characteristics of a sandy loam soil, since the soil texture analysis had defined the majority of the soil samples at different depths as sandy loam. The average divergence of the predicted SWR curve for the A p horizon (Fig. 9c) from laboratory measurements was -13%. The largest differentiation occurred  at -1000 cm H 2 O (Fig. 9), while for near-saturation (-10 cm H 2 O) and -15 cm H 2 O the deviation was -1.5 and -0.9%, respectively. However, in the model simulations, the matric potential only decreased to values near -1000 cm H 2 O during the dry periods in the calibration. Therefore, these discrepancies from the measured data were assumed to not significantly affect the results. The SWR curves for the other three soil horizons similarly followed a typical curve of a sandy loam soil (Karup et al., 2017;Pittaki-Chrysodonta et al., 2018).
The predicted soil water characteristics for all the soil horizons were also evaluated by the Campbell b parameter, which is a pore-size distribution index equal to the negative slope of the soil-water retention curve on a log-log scale [log(-h) vs. log(q) system] (Campbell, 1974). Typical values of Campbell b for coarse-textured sandy soils (macropores) are in the range of 1.0 to 5.0, while for medium-textured soils (sandy loam) the range is 5.0 to 14.0 . For this study, the Campbell b parameter for the macropores varied from 1.4 to 1.6, while the calculated Campbell b using the estimated SWR curve parameters from the single-porosity model yielded to 6.0 to 13.2 and therefore these ranges are consistent with the literature review.
Water levels measured by the piezometers during the simulation periods for all models were compared with the groundwater level measurements from P18 (Fig. 10). The P18 piezometer was selected as the most representative well for the drainage subcatchment, since P17 and P19 were located at the border of the field outside the tile-drained area (Fig. 1). The measurements of the groundwater level from P18 had moreover shown poor correlation with the measurements from P17 and P19. For the calibration model, the dual-permeability model had shown strong correlation between the measured against the predicted values (r = 0.98) and thus a good ability to describe the water level dynamics pattern of the specific agricultural landscape. In contrast, the dual-porosity model had a poorer correlation (r = 0.88) and consequently resulted in a predicted pattern that was less close to observations. The lower predictive ability related to the water level dynamics behavior of the dualporosity model was expected due to the restriction of water flow to macropores and no flow between aggregates. Therefore, the lower predictive ability was also related to the flatter shape of the SWR curve and the high K s,f and the instant transitions from saturated to drained conditions. Moreover, the assumption of the similarity in soil water retention curves between the two pore regions limited the capability of the immobile domain to sustain water. These characteristics were the reason for the abrupt termination in the recession limbs of drainage for the dualporosity model, revealing the change to unsaturated conditions at the depth of the tile drainage. The single-porosity model, on the other hand, showed a good ability (r = 0.94) to simulate the observed drainage dynamics pattern. The first validation model showed similar results to the calibration models (Fig. 10). For the second validation model no comparisons were made between the simulated and observed groundwater level due to the shorter time of simulation and, consequently, smaller number of observed measurements. Despite the fact that the dual-permeability and single-porosity models had shown strong correlation to the measured groundwater data, the slope of the regression analysis was close to 0.64, indicating an underestimation of approximately 36% from the 1:1 line (Fig. 10). However, P18 is located close to the depression zone and observed values of the groundwater level were therefore relatively high compared to the areas outside the depression zone. Moreover, the estimated water level from the three approaches can be regarded as average water levels for the entire drainage subcatchment. Due to the hydro-topographical gradients, the non-systematically drained field and the spatial variability in the soil hydraulic properties, the observed local water levels at different points in the field are expected to vary compared to the predicted values. Nevertheless, the water level output from the simulations with the dual-porosity and especially the dualpermeability model reproduced well the behavior of the water level dynamics occurring in the depression zones based on the calculated correlation coefficients. This is a strong indicator of the decisive influence of flow and transport pathways in depression zones on the soil-water balance and drainage dynamics in the studied area.

Comparison of the Three Approaches
For the calibration as well as for the second validation dataset, the dual-permeability model was better able to describe the observed drainage dynamics in the drainage subcatchment, while for the first validation dataset the best results were obtained with the dual-porosity model. All the models failed to achieve a good fit during the dry periods, probably due to an unknown source (e.g., baseflow) that was recharging the tile drainage network. However, the failure to predict the low summer drainage discharge has no quantitative impact in terms of estimating nitrate transport by drainage (Neal, 1991). Furthermore, this study aimed to examine the capabilities and the weaknesses of the selected models. Depending on the processes that are the target for the modeling, the appropriate modeling approach can be selected.
The simultaneous, independent calibration of 28, 17, and 4 parameters of the dual-porosity, dual-permeability and singleporosity models, respectively, could lead to non-unique solution sets. Moreover, the initial parameterization could have influ-enced the final estimation of the effective parameters. However, the validation of these parameters with data obtained by the laboratory or in situ measurements of the groundwater level showed that the single-porosity and especially the dual-permeability models were capable of representing satisfactorily the average actual internal flow and transport pathways within the studied integrated hydrogeological system. This result may be useful for other studies involving the modeling of nutrient transport where an accurate representation of the physical processes in the flow domain would be important to simulate the residence time of nitrate, the process of nitrate reduction in the subsurface and the evolution of the contaminant transport.

CONCLUSIONS
The aim of this study was to examine the ability of different model approaches to simulate the water level and tile drainage dynamics in an agricultural field characterized by intense hydrotopographical gradients and drainage sequences in a marine catchment. Three models in HYDRUS-2D were used having as a common feature the inclusion of preferential water flow. Despite the differences between the conceptual approaches, all models produced a reasonably good fit to the hydrograph features for both the calibration and validation datasets. The dualpermeability model was best able to predict satisfactorily the drainage dynamics in the selected landscape. However, it is complex and has a large number of parameters needing calibration, whereas the single-porosity model gave an acceptable accuracy with significantly fewer parameters to calibrate. There is a need for further study to determine the exact type and magnitude of the unknown source that recharges continuously the installed tile drainage network and produces small discharge rates during the dry periods. There is also a need for further studies to examine the accuracy of the calibrated saturated hydraulic conductivity of the macropores and, consequently, the indirect estimate of the potential water flow through the preferential flow pathways. To this end, in situ measurements of soil water content would be a significant aid to achieve a more reliable, calibrated water flow model. An alternative modeling approach could be a three-dimensional, dual-permeability model to take into account the spatial variability in the geological layering across the entire agricultural landscape, the lateral flow velocities in macropore domain, and the characteristic of non-systematically, tile-drained fields. A tracer test using a conservative tracer and simulation of solute transport would facilitate closer examination of the transport pathways.
Finally, an integrated, three-dimensional, hydrobiogeochemical model could be developed to simulate the solute transport and the nitrogen dynamics in the subcatchment using as input the outcome of this study, which will be the average soil water release characteristics and the hydraulic conductivity for the two different transport pathways.