Important factors when simulating the water and nitrogen balance in a tile-drained agricultural ﬁ eld under long-term monitoring

• Daisy model simulates tile drainage discharge and nitrate leaching under long term crop rotation satisfactorily • Freezing and thawing processes in the soil cause the major discrepancies be-tween simulated and measured nitrate leaching • Carry-overeffectofthewinterrapecrop leads to higher than measured nitrate leaching to drains by the model • Highnitrateconcentrationsinthedrain-age water signify the need for better management practices


Introduction
Nitrogen (N) is one of the primary and essential plant nutrients.Through its repeated cycle, which consists of non-sequential steps of fixation, assimilation, ammonification, nitrification, and denitrification, N circulates in various form in nature.N plays a vital role in agricultural systems for achieving their maximum economic yields and high quality products.However, the increased application of N fertilisers has had several consequences for the environment and raised concerns about the effectiveness of N inputs (Fowler et al., 2013;Baram et al., 2016), resulting in the agricultural sector being the main contributor to the nonpoint pollution of surface and groundwater resources (Bakhsh et al., 2000;Bakhsh et al., 2001;Bakhsh et al., 2004;Beaudoin et al., 2005;Billy et al., 2011;Mollerup et al., 2014).Nitrate (NO 3 − ) is the main source of aquifer contamination in agricultural areas (Lee et al., 2020), where the use of excessive amounts of N can result in a high degree of NO 3 − being transported via water pathways, causing accelerated eutrophication of lakes and rivers (Erisman et al., 2013;Dupas et al., 2015).
The presence and impact of tile drain systems can play a decisive role in the quality of aquatic ecosystems, with tile drains facilitating the rapid transport of nutrients and dissolved chemicals, and hence bypassing the soil domain (Jaynes et al., 2001;Dinnes et al., 2002;Van den Eertwegh et al., 2006;Lapen et al., 2008aLapen et al., , 2008b;;Syswerda et al., 2012;Cordeiro and Sri Ranjan, 2015;Amado et al., 2017;Schilling et al., 2019).However, tile drain systems are crucial in naturally poorly drained agricultural soils as important water management practice that provide efficient crop production systems (Tiemeyer et al., 2010;Nangia et al., 2010).By controlling the groundwater table and removing excess water in the root zone, these systems enhance crop growth (Sloan et al., 2016) and prevent surface runoff.
In Denmark, more than half of the country's agricultural fields are artificially drained (Olesen, 2009;Møller et al., 2018) and thus accurate estimates of drainage and nutrient transport are of great importance for policy support and the development of improved crop management systems.Several factors affect N leaching from the surrounding soil to the tile drain system, such as precipitation, temperature, soil characteristics, the presence of macropores, crop varieties and field management (Ernsten et al., 2015).According to the Nitrates Directive (91/676/EEC; EEC, 1991) and the Drinking Water Directive (98/83/EC), the European limit for NO 3 − -N concentrations in drinking water supplies and coastal water is 11.3 mg N L −1 .To minimise N losses to water produced for drinking and hence be able to ensure NO 3 − -N concentrations below that limit, more effective N management strategies are needed that are suited to local and regional characteristics.Determination of N fate is still a challenge and requires comprehensive field operations and measurements, which are costly and labour intensive.Furthermore, the data obtained from monitoring programmes rarely reveal the presence of dominating fate processes and therefore do not explain the entire N balance.To address this, modelling tools that incorporate important N-fate processes and are capable of capturing water and N balances play a crucial role in understanding these complex natural systems, developing management strategies and subsequently mitigating N losses to the environment.
In the last few decades, there has been a focus on crop rotation modelling as an effective tool for investigating nutrient losses over extended periods among certain crop sequences (e.g.Kovács et al., 1995;Blombäck et al., 2003;Post et al., 2007;Beaudoin et al., 2008;Salado-Navarro and Sinclair, 2009;Syswerda et al., 2012;Hlavinka et al., 2014;Sun et al., 2017;Gillette et al., 2018) and for studying the development of biomass, crop dry matter and N yields (Berntsen et al., 2006;Nendel et al., 2013;Negm et al., 2017).However, only a few studies have included N leaching to artificial drains (Conrad and Fohrer, 2009;Marjerison et al., 2016;Negm et al., 2017;Liang et al., 2018;Gillette et al., 2018;Matteau et al., 2019;Yin et al., 2020).Using the CoupModel, a physically based soil-vegetation-atmosphere transfer (SVAT) model, Conrad and Fohrer (2009) calculated NO 3 − -N leaching from an artificial subsurface drainage system using data from an eight-year monitoring programme in a North German field under winter wheat/red clover rotation and organic farming conditions.The CoupModel showed discrepancies in drainage and NO 3 − -N leaching to the tile drain system for short-term events during the seepage period, and an intensive parameter adjustment was necessary for the second period with red clover crop (Conrad and Fohrer, 2009).Matteau et al. (2019) coupled an empirical NO 3 − production model adapted to the regional conditions of a study area near Québec City with the one-dimensional HYDRUS model to predict NO 3

−
-N leaching into subsurface drains under maize, canola and spring wheat crop rotation for 2009 and 2011.Although the simulation results for the drainage were stated to be satisfactory, the lower boundary of the model matrix was considered free drainage, with the assumption that all water was recovered by the subsurface drain (Matteau et al., 2019).Therefore, the groundwater dynamics and NO 3 − leaching below the drain depth were not investigated.Additionally, this model study did not offer any insight into the long-term effects of fertiliser inputs on N leaching as it only applied data from 2009 and 2011.In contrast, Yin et al. (2020) simulated N leaching using the STICS model in a 34-year crop rotation experiment with winter wheat, sugar beet, oil seed rape and spring barley in rotation.In that study, the model proved its ability to simulate aboveground biomass by capturing the measured crop N uptake but was not able to capture the high N leaching events to the drains, reported to be due to limitations in calculating the fate of accumulated N in deep dead roots (Yin et al., 2020).Additionally, Marjerison et al. (2016), who used the PNM model to simulate drainage and NO 3 − leaching from an artificially drained maize field in two fields in New York and Minnesota, USA from 1994 to 1999, reported discrepancies between the simulated and measured values due to unaccounted preferential flow and the PNM model's inability to cope with high precipitation events or snowmelt in spring.Liang et al. (2018) also reported discrepancies in simulated daily drainage discharge by the process-based water and N management model WHCNS due to unaccounted preferential flow, and ignored microclimate conditions, while Kollas et al. (2015) compared the ability of fifteen crop-growth simulation models to predict crop yields in rotation at five sites across Europe and described the high quality experimental data as the "backbone" for the development and improvement of multi-year simulation models.Furthermore, the study of Kollas et al. (2015) highlighted the three challenges that needed to be addressed in future studies: i) the carry-over effect of crop rotations on nutrient balances in agro-ecosystem modelling approaches, ii) the inclusion of adequately investigated crop properties (e.g. for rapeseed oil and sugar beet) in crop rotation datasets and models, and iii) the identification of frost effects on intermediate crops in long-term modelling studies.Thus, in the present study a one-dimensional Daisy model (Hansen et al., 1991) was utilised to simulate the water and N balance while addressing under-investigated issues in modelling such as the effects of snowmelt and frozen soil conditions on water and solute transport, the carry-over effect of fertilisers and crop residues on NO 3 − -N leaching, the dynamics of NO 3 − -N leaching to drains and groundwater, the groundwater table dynamics at field-scale with tile drain systems and, finally, the degree of preferential flow and fast transport of nutrients through the soil.The Daisy model has been used in a number of studies to predict the water balance and/or NO 3 − leaching at catchment scale (Styczen and Storm, 1993;Refsgaard et al., 1999;Thorsen et al., 2001;Boegh et al., 2004;Hansen et al., 2007;Van der Keur et al., 2008;Hansen et al., 2009).However, to the authors' knowledge, rather than the PhD thesis of Nagy (2020) focused on the degree of preferential transport, no studies have simultaneously investigated N uptake by the crops and NO 3 − leaching out of the root zone to the tile drain system and groundwater using the Daisy model on a long term crop rotation dataset.Moreover, none of the abovementioned components of the N balance has been investigated in combination with the importance of preferential flow and realistic estimates of N losses through denitrification processes.Therefore, using comprehensive monitoring data over a ten-year period from a tile-drained agricultural loamy soil under a longterm crop rotation, the aims of this model study were: i) to determine the dominating processes that are essential for N-leaching risk assessments, ii) to estimate the water and N balance of the field, including the quantification of components that are difficult to monitor/measure, and finally iii) to investigate the causes of potential discrepancies between measured and simulated water and N balance data with a focus on NO 3

−
-N leaching to a field's tile drain system.

Field and crop management description
The tile-drained agricultural field Faardrup (2.33 ha; 55°19′ N, 11°2 0′ E, elevation 28 m above sea level) is located south of Slagelse on Zealand, Denmark, and has since 1998 been included in the Danish Pesticide Leaching Assessment Programme (PLAP) (n.d), Rosenbom et al. (2015) http://pesticidvarsling.dk/?lang=en.It consists of sandy loam (Argiudoll mapped with both Hapludoll and Vermudoll, based on USDA Soil Taxonomy) and is situated on a clay till plain from the Late Weichselian age with a small slope (0-3%).
Subsurface tile drains were installed in the field between the 1940s and 1960s.The drains were established at a depth of approximately 1.2 m and with a horizontal spacing of approximately 18 m.Before monitoring began, the field's tile drain system was isolated from the surrounding tile drain systems to ensure that drainage only came from the 2.33 ha field.Additional information about the site characterisation and experimental design is described in Lindhardt et al. (2001).The farming practice in the field is in line with the conventional practice of the region and with applications of N as recommended as good management practice in the Danish Action Plans for the Aquatic Environment during that period (Ministry of Environment and Food, 2017).Each year, the amount of N fertiliser applied has been adjusted to suit the selected crop type and previous year's climatic conditions (Ernsten et al., 2015).Spring barley and winter wheat were the most commonly grown crops during the monitoring period 2000-2009 (Table 1; Ernsten et al., 2015;Rosenbom et al., 2015).Detailed data about crop cultivation and development (sowing, emergence and harvest time), fertiliser input and field operations were used to create the management file as an input to the Daisy model.The sowing and harvest dates of each crop in the rotation and the amounts and dates of N fertiliser (inorganic only) applied from 2000 to 2009 are presented in Table 1.

Instrumentation and measurements
Automated water samplers (Teledyne ISCO Lincoln, NE, USA) were used to collect drainage samples and the drainage water was sampled time proportionally until July 2004 and flow proportionally for the rest of the period.In relation to the flow proportional sampling, subsamples of 200 mL were collected for every 3000 L of drainage water during the winter season (September-May) and for every 1500 L during the summer season (June-August).The collected subsamples were pooled, and a subsample of these transferred to the laboratory on a weekly basis and refrigerated until analysis.NO 3 − concentrations were measured in each drainage sample by applying the method described in Ernsten et al., 2015.The total mass or flux of NO 3 − -N (kg N ha −1 ) transported to the outlet of the tile drain system was calculated by multiplying the measured concentration of NO 3 − -N by the drainage volume monitored since the previous collected drainage sample, considering the size of the treated area of the field.The water table at Faardrup was measured every other week with a hand-held water-level meter in four piezometers installed in each corner of the field (Fig. 1).Two groups of timedomain reflectometry (TDR) probes were installed in 1999 in the northwest corner of the field for hourly measurement of the volumetric soil water content at depths of 0.25 m, 0.6 m, 0.9 m, 1.1 m and 2.1 m.At each depth, three probes were installed and the average of these is applied in this study.Furthermore, soil water from the variably saturated zone was extracted for analysis of NO 3 − -N concentration by applying the method described in Ernsten et al., 2015.Here, in 1998 eight suction cups were installed in clusters on each side of the field about 20 m from the northwest corner.In each cluster, four suction cups were installed at depths of 1 m and 2 m (Lindhardt et al., 2001;Ernsten et al., 2015).

Daisy model
Daisy is a physically based model developed by Hansen et al. (1991) for mechanistic simulations of the water and N balance of agricultural fields.The main modules and processes described in Daisy include water flow, plant water and N uptake, soil heat and N turnover and transport dynamics.Water flow in the soil matrix is described by Richards' equation and solute transport is described by the convection-dispersion equation.Steady-state tile drainage flux per unit surface area is computed by Hooghoudt's (1940) equation, allowing only water to flow towards the tile drain when the groundwater table is above the tile drain (Mollerup et al., 2014).Macropores can be introduced as a fast flow and transport domain, which in Daisy is defined as the transport of water and solutes in large biopores or fractures where the capillary forces no longer dominate.The macropore module in Daisy, designed and described by Mollerup (2010), has been applied as the preferential flow domain in both the 1D and 2D versions of Daisy (Abrahamsen, 2011), and has been tested in technical reports published by the Danish Environmental Protection Agency (Hansen et al., 2010a(Hansen et al., , 2010b)).The macropore model divides the fast flow domain into a number of classes, instantly transporting the water either to the bottom of the macropore, where it might slowly leak back to the soil matrix domain (matrix macropore), or to the tile drain (drain-connected macropore).The main properties applied by the macropore module are the start depth/point, end depth/point, density and diameter of the macropore.For the drain-connected macropores, the position or installation depth of the tile drain must be specified.For the matrix class macropores, the hydraulic conductivity of the macropore wall is specified relative to the hydraulic conductivity of the surrounding soil matrix.The wall thickness of the macropores is assumed to be 10% of the diameter of the macropore (Abrahamsen, 2011;Abrahamsen, 2020).Moreover, the density of the macropores is defined as a function of the distance to the tile drains (Abrahamsen, 2011;Nagy, 2020).Soil organic matter turnover is based on multiple pools with different turnover rates, which are affected by soil heat, moisture and clay content (Hansen et al., 1991;Hansen et al., 2012).In the crop model, dry matter is partitioned into storage organ, leaves, stem and roots.Photosynthesis is linked to a soil-vegetation-atmosphere transfer (SVAT) model of water and energy, which is affected by water and nitrogen stresses.Crop water uptake is calculated by the Darcy flow based on climatic demands and root crown water potential.Nutrient uptake of the crop is based on both convection and diffusion (Hansen et al., 1991;Abrahamsen and Hansen, 2000;Hansen et al., 2012).In the present study, the 1D Daisy model setup was used to simulate drainage discharge (Q) with and without the macropore model, groundwater table (GWT), NO 3 − -N leaching (N.Leach) through the matrix (with and without the macropores) to the tile drain, soil water content (SWC), crop dry matter (DM) and crop N (CrN) yield.Moreover, denitrification rates for each season were simulated as a component of the total N balance.

Meteorological data
The meteorological data used in the current study were recorded at the meteorological station at Aarhus University Flakkebjerg, 4 km east of the Faardrup field.As an input to the Daisy model, a weather file was generated where measured precipitation (P), average air temperature (T), global radiation, relative humidity and wind speed were provided at a daily time step.The model calculated the potential evapotranspiration (ET) using the Penman-Monteith equation (Monteith, 1965).Using the correction factors of Allerup and Madsen (1980), P was corrected to the ground truth or surface values.In addition to these meteorological data, nine years of meteorological data from Flakkebjerg covering the period 1991-1999 were used as a warm-up period for the model, including the actual registered crop history for the same period, to minimise the effects of the initial values of organic matter pools on the results (Styczen et al., 2005).

Hydrogeological setting
Saturated hydraulic conductivity (K s ) and soil water characteristics at different horizons were determined based on the pedological description and extensive measurements on large soil cores (Iversen et al., 2004) from three two-metre deep soil profiles in the northern, eastern and western sides of the field.Iversen et al. (2004) used 100 cm 3 and 6280 cm 3 soil core samples to determine the soil hydraulic characteristics based on the VG-M model (n, α; Table 2) (Van Genuchten, 1980) and hydraulic conductivity (K s ) respectively.An overview of the physical properties of the three distinct soil horizons (A, B, and C) represented in the soil profiles is presented in Table 2.As part of the PLAP project, a MACRO model (Larsbo and Jarvis, 2003) was set up and calibrated based on the above-mentioned estimated hydraulic parameters, groundwater table depth, soil water content at three depths, drainage, bromide concentrations in suction cups, and bromide concentrations in drains for the period May 1999-May 2004.The model results were validated against water balance and bromide transport for a two-year period (Barlebo et al., 1999;Kjaer et al., 2007;Rosenbom et al., 2015).Based on these calibrated hydraulic parameters, and the lithological descriptions provided by Lindhardt et al. (2001) and Iversen et al. (2004), a five-metre deep soil profile n, was set up in the Daisy model.Soil properties from 185 to 500 cm were considered as homogenous and equal to the C horizon.In a study of long-term pesticide leaching at Faardrup, Rosenbom et al. (2015) indicated that the presence of preferential transport through macropores is likely (burrows and fractures) in this field.Therefore, based on the excavations and profile descriptions of Lindhardt et al. (2001), both single and dual permeability scenarios were set up with the Daisy model and tested.These are described in section 3.3.2below.

Lower boundary conditions
Hooghoudt's drainage theory was initially designed for lateral water flow into tile drains under an impermeable lower boundary condition (Hansen et al., 1991;Mollerup et al., 2014).The Daisy model enables the user to choose other lower boundary conditions, such as an aquitard or an actual measured groundwater table, provided as an input file within the temporal resolution of the weather input data.Defining an aquitard as the lower boundary condition results in a much faster model computationally (Mollerup et al., 2014).The two main parameters of the aquitard lower boundary condition in Daisy are the hydraulic conductivity of the aquitard (k_aq) and the thickness of the aquitard (z_aq).Thus, an aquitard with default hydraulic parameters based on the Daisy model manual (Abrahamsen, 2020;Mollerup et al., 2014; k_aq = 0.001 cm day −1 and z_aq = 200 cm) and a constant soil water pressure of −500 cm H 2 O at the lower boundary, corresponding to the bottom of the model domain, were used in the current study.

Measures of accuracy and variable correlation
The accuracy of all the time-series simulations (Q, SWC, GWT, and N. Leach) was measured by the metrics of the Nash-Sutcliffe efficiency (NSE) and the Kling-Gupta efficiency (KGE): where X simi is the predicted value for the i-th date, X obsi is the observed or measured value for the i-th date, X obs is the mean of observed or measured values, and n is the total number of days throughout the simulation period, and.
where σ obs is the standard deviation in observations, σ sim is the standard deviation in simulations, μ sim is the simulation mean, and μ obs is the observation mean (i.e.equivalent to X obs ).
Based on Knoben et al. (2019), NSE and KGE values should not be directly compared as their relationship is non-unique and to a high extent depends on the coefficient of variation of the measured time series.The KGE measure of Gupta et al. (2009) seems able to overcome some of the limitations of NSE, such as high and low flow values that are not accurately represented (Pechlivanidis et al., 2010).Mean flow benchmark, which is an inherent feature of NSE, results in NSE equal to zero and, for most modellers, this is the threshold for the performance of the model.However, KGE does not have an inherent threshold and the mean flow benchmark equals −0.41 for KGE.Therefore, values between −0.41 and 1 could be considered a "reasonable" model performance (Knoben et al., 2019).For both criteria, 1 indicates a "perfect" match between the observed and simulated values (Nash and Sutcliffe, 1970;Gupta and Kling, 2011), therefore modellers must be specific about the benchmark against which they are comparing their model performance (Knoben et al., 2019).
To measure the robustness of the model for the simulations of DM and CrN, the root mean squared error (RMSE) was used: To assess the dependence of various parameters within the model and to determine the strength of a relationship between simulations and measurements, a Spearman's rank-order correlation test was carried out.This method has been used in multiple studies to understand the interrelation of model parameters and the dependence/independence of the data in ecosystem and agricultural modelling (Callesen et al., 2007;Härdtle et al., 2007;Weymann et al., 2008;Wójcik-Gront, 2020).Spearman's rank correlation coefficients have values between −1 and 1.Values closer to −1 indicate a strong negative correlation and values closer to 1 mean a strong positive correlation between two parameters.

Results and discussion
This section presents the simulation results on the water and N balance compared with the actual measurements and discusses the model performance and discrepancies.

Water balance
Measured daily drainage dynamics (Q) from 2000 to 2009 showed that the major events during the period occurred in the winter half-year (Fig. 2a), with the highest daily event (15 mm) in February 2007, while during the summer almost no flow was measured.A high evapotranspiration due to an increased plant water update in the top 1 m of the soil during summer are the main reasons for significantly lower drainageQ during this period.Simulation of daily Q throughout the 10-year period (Fig. 2a) resulted in NSE and KGE values of 0.27 and 0.61 respectively.In view of NSE being highly sensitive to extreme values, the advice being to account for any outliers (Moriasi et al., 2007) and KGE being within the satisfactory range (Fig. 2a), simulations were analysed for each year separately to investigate the reasons behind the low value of NSE for the entire 10-year simulation period.Based on the simulation results for the cumulative Q (Table 3), the model performed satisfactorily throughout the simulation period except during the hydrological year (from 1 October to 30 September each year) 2003-2004, when the model overestimated Q by 42 mm (Table 3).The model's results of the cumulative Q for the entire simulation period (Fig. 2b) demonstrated a perfect match between simulated and measured values based on NSE a satisfactory match based on KGE.Based on Moriasi et al. (2015) and Muma et al. (2016), the results could be classified as highly accurate and satisfactory for field-scale models, especially when simulating tile drain discharge including summer/dry periods when Q events rarely occur.Furthermore, compared with the results reported for simulations of drainage discharge by the HYDRUS model (Matteau et al., 2019) and the CoupModel (Conrad and Fohrer, 2009), NSE obtained for the 10-year simulated Q with the Daisy setup model was higher.
In order to analyse water dynamics in the vadose zone, simulations of SWC by the model were compared with TDR measurements with probes installed at various depths (Fig. 3a, b, c, and d).The simulated SWC at 25, 60, 90, and 110 cm depth showed good agreements overall with the TDR measurements.However, the SWC for all depths (Fig. 3a, b, c, d) was underestimated by the model in 2007 and 2008, and overestimated for the shallower depths in 2003-2004 (Fig. 3a and b).NSE tended to increase with depth, reaching the highest value of 0.60 at the depth of 90 cm (Fig. 3c) and the lowest value at the depth of 60 cm (Fig. 3b).Conversely, the highest KGE (0.73) was obtained for the simulations of SWC at the depth of 90 cm (Fig. 3c) and the lowest KGE (below the threshold) at the depth of 25 cm (Fig. 3a).In general, the one-dimensional Daisy model with one set of defined hydraulic parameters for each horizon demonstrated a good performance in simulating the SWC in the top 1 m of an expected highly variable three-dimensional tile drained field where in reality the hydraulic parameters for each horizon could differ quite markedly from the calibrated parameter set.
Simulated GWT (Fig. 4) was compared both with the average measured values obtained from the four piezometers located in each corner of the field (P1-P4, Fig. 1) and the measurement from P1, which is the piezometer closest to the drainage outlet.Values of NSE and KGE (calculated for the 15-day resolution based on the measurements) were very similar in both comparisons for the 10-year period (Fig. 4).The model simulations, however, showed a slightly better match with the average values of P1-P4 (Fig. 4).Measurements from P1 showed a higher winter water table at the start of the simulations (2000 to 2002) as also described by the model.However, in 2007 and 2008, simulated GWT during the summer was much lower than the measurements (Fig. 4).Hooghoudt's drainage function is based on the assumption that the GWT is simulated between two tile drains, and the measurements  (Van Genuchten, 1980) for the three horizons included in the Daisy model (Lindhardt et al., 2001) collected in the buffer zone did not reflect these exact circumstances.Optionally, a second soil column could be parameterized in the model for the buffer zone specifically, where the distance between drain pipes and the position of the calculation node are adjusted, to assess the dynamics of the water table.However, we aimed at representing the dynamics within the entire field and the abovementioned optional step would have not contribute to the better understanding of the system.Therefore, it was concluded that the GWT simulations by the Daisy model were able to adequately reproduce the measured data over the entire period.The agreement between both the average and P1 measurements with the GWT simulations, with NSE of 0.63 and 0.60 respectively, showed a better performance by this Daisy implementation than the models used in similar studies, such as the CoupModel (Conrad et al., 2009).An overview of the average annual water balance simulated by the model from 2000 to 2009 is presented in Fig. 5.The annual water balance for each year is calculated on 30 September, the day which the hydrological year ends, and finally averaged for the entire 10-year period.The average annual P input to the model was 675 mm (hydrological year from 1 October to 30 September), while the annual average output consisted of 466 mm ET, 118 mm deep percolation, 97 mm Q, and 6 mm soil water storage change.Storage describes the differences in water content in the soil profile from one year to another year.For example, this could be related to a wet or dry spring or differences in groundwater levels in each year.

N balance
NO 3 − -N concentration measurements in the soil water at 1 and 2 m depth and the drainage at 1.2 m depth collected from the field are presented in Fig. 6.Except for some periods from 2007 to 2009, concentrations measured in the water collected from the suction cups at both depths (Fig. 6a) were generally above the 11.3 mg L −1 EC limit (European limit for NO 3 − -N concentration in drinking water and coastal water defined by the Drinking Water Directive; 98/83/EC).In agreement with the suction cup measurements, measured NO 3 − -N concentrations in the drainage samples from 2000 to 2009 (Fig. 6b) were also mostly above the EC limit.The highest concentrations were measured from 2000 to 2002, when winter wheat and sugar beet crops respectively were grown, and from 2005 to 2006, when maize and spring barley respectively were grown.For both the above-mentioned crop rotations, there was a seven-month period between the two crops where the soil was left bare.The NO 3 − -N concentrations measured within these periods were very high.Ernsten et al. (2015) also indicate that the highest NO 3 − -N was measured after crops with a short growing season followed by bare soil.In agricultural systems with rotating annual crops, leaving the soil bare for a period or part of the year will lead to an increased leaching of soil N due to the lack of uptake by the plants (Sparks, 2014).The source of this N is usually the organic residuals of the previous crop season, which are then mineralized and leached with the water in the matrix through the soil column.In the period from 2000 to 2009, the Danish government imposed regulations and restrictions on the use of chemical   fertilisers in agriculture, therefore the inputs during this study (Table 1) were sub-optimal (Dalgaard et al., 2014).Nonetheless, NO 3 − -N concentrations in the drainage indicated the need for improved management practices and strategies to decrease leaching to a level below the EC limit.
Further studies should be carried out to reveal the trend in N losses under current farming practices.
Daily N.Leach simulation results showed an acceptable match compared with the measured values (Fig. 7a).However, the results showed  lower values of NSE and KGE (0.27 and 0.42 respectively) compared with the cumulative leaching.Cumulative N.Leach via the tile drain system over ten years was simulated satisfactorily with NSE of 0.87 and KGE of 0.74 (Fig. 7b).The main divergence in the model was in 2004 when the simulated leaching was approximately 20 kg higher than that measured.Furthermore, in 2007, simulations showed a 10 kg higher N.Leach compared to the measurements, totalling 30 kg from 2000 to 2009 (Table 4; Fig. 7).Simulated CrN and DM for each season in the 10-year rotation period were compared with the measured values (Fig. 8a and b).The results indicated a good match, with RMSE for CrN and DM of 14.12 kg N ha −1 and 1.23 Mg DM ha −1 respectively.Accurate model simulations of crop N uptake (CrN in Fig. 8a) are extremely important for simulations of the entire N balance, especially in relation to correct estimation of N leaching to tile drains and deeper through the soil profile towards the groundwater.Due to some technical issues during the harvest of sugar beet in 2009, DM was not measured and the applied value was calculated based on the recorded fresh weight of the harvested roots and the common standards for dry matter percentages in sugar beet.Based on Plauborg et al. (2010), the accuracy of DM simulations compared with actual measurements is an  indirect indication of the model's realistic estimations of ET.Evidently, the overall fit between the simulated and measured DM could be categorised as satisfactory (Fig. 8b).
In view of the fact that proper estimates of N transformation processes provide the basis for modelling N losses (Thorp et al., 2006;Qi et al., 2011;Negm et al., 2014;Liang et al., 2018), an annual average was calculated for different components of the N balance based on both measured and simulated values in the period from 2000 to 2009 (Fig. 9).Based on the simulations, a total amount of 143 kg N ha −1 was denitrified during the period, corresponding to an annual average of 7% or 14 kg N ha −1 .Furthermore, 71% (139 kg N ha −1 ) of the total average annual N input was utilised by the crops and harvested at the end of the season, while 11% (22 kg N ha −1 ) was lost due to deep percolation and 9% (17 kg N ha −1 ) was leached through the tile drains to the surface water.A yearly summary of the N budget and main processes simulated by the model for the crop rotation in 2000-2009 is presented in Table A.1.

Issues related to precipitation and crop parameterisation
In 2004 (Fig. 10b), two Q events were simulated in January-February and March-April that were not recorded in the measurements to the same extent.The events led to a difference in Q of approximately 50 mm in the spring period and a total yearly difference of 42 mm (Fig. 10b and Table 2 respectively).There could be several reasons for the discrepancy between the simulated and measured values.One could be that the P was measured at a meteorological station four kilometres away from the field, and thus it did not fully represent the P at the field, especially during localised heavy rain events.Another reason could be the model's inability to simulate slower infiltration rates of P or snowmelt under frozen soil conditions.Below-zero temperatures in January meant that the upper soil water froze and probably remained that way into February, despite air temperatures being above zero in the second half of January.As presented in Fig. 10a, the model simulated a considerable amount of snow storage following the below-zero temperatures.During this period, the soil would have had a lower hydraulic conductivity as the soil water was frozen, whereas the model calculated soil water infiltration and flow with unchanged and higher hydraulic conductivities (Ireson et al., 2013;Holten et al., 2018).Furthermore, the model peak drainage flows not observed in the measurements can be explained by the modelled infiltration of snowmelt water, which could not take place in reality due the frozen soil conditions.Hence, the model did not capture the natural low infiltration of water through the soil from January to April.Instead high and sudden unmeasured events were simulated in January and February.With increasing T above 0 °C at the start of April 2004, a high event of Q was measured  in connection with a short period of P.This could most likely be explained by the increasing temperatures and subsequent thawing process that started due to temperatures above 0 °C.Hence, frozen soil water was suddenly released and initiated the downward flow, which increased with the simultaneous occurrence of 6 mm of rainfall.The high event in measured Q and difference from the model simulation could be explained by the fact that soil water processes such as freezing and thawing were not included in the model.Considering the model's inability to cope with frozen soil conditions and, consequently, the lowered hydraulic conductivity during the frost period, this hypothesis offers the most likely explanation for the discrepancy between measured and simulated drainage events from January to -April 2004.However, the total simulated and measured accumulated drainage from January to May would be expected to be the same, which strongly suggests that a higher value of P was measured in that period at the meteorological station at Flakkebjerg than at the field site.Consequently, N leaching to the tile drain system accumulated throughout the year was simulated to be approximately 20 kg N ha −1 higher than measured (Fig. 10c).Based on this overestimated N loss by the model in the beginning of 2004, a simulated N stress in the crop (winter wheat sown in September 2003 and harvest in August 2004) was expected.However, the simulation results did not indicate any days of N stress during the winter wheat growing season.Therefore, NO 3 − -N concentration measurements from the suction cups at 1 and 2 m depth were compared with concentrations simulated by the model at these depths (Fig. 11a and b).The model showed a significantly higher NO 3 − -N concentration in the soil at 1 m, which was not measured (Fig. 11a).Further analysis of the Daisy output indicated that the source RMSE=1.23 RMSE=14.12  of this simulated accumulation of N at this depth was most probably from the crop residuals in the soil after the harvest of winter rape in July 2003, which was the same amount of 20 kg N ha −1 .In the model, mineralized residuals converted to NO 3 − -N were not taken up by the following crop.The roots of the simulated winter wheat in 2004 penetrated to a much lower depth than the same crop in 2000 (80 and 120 cm respectively) as they were limited by the high groundwater table during the winter, leaving the accumulated N out of the roots' reach.Therefore, due to the incorrect P input to the model and/or the model's inability to handle frost and thaw processes, it overestimated the N concentration at 1 m depth and hence resulted in an overestimation of N loss to the tile drain system.Finally, it is worth noting that the sampling of drainage for NO  In late February 2007 (Fig. A.1b), the model was unable to capture a high Q event, which again could be related to the use of P from the climate station situated 4 km away from the field and/or the model's incapability of coping with soil frost and lowered infiltration rates (Fig. A.1a).Consequently, the simulated N.leach from January 2007 to March 2007 was lower than that measured (Fig. A.1c).Based on the findings of Holten et al. (2018), who investigated the effect of freezing and thawing on drainage and solute transport in loam and silt soils, significantly higher volumes of water percolate from frozen topsoil in loam and for both soil types, and with the start of thawing in the frozen soil columns, the ponded water infiltrates quickly and percolates through the soil.This effect could be potentially responsible for the underestimation of Q by Daisy (Fig.

Water and solute transport through macropores
The degree of preferential flow and the input of the macropores to Q and N.Leach transport were investigated in a dual-permeability scenario (AM: active macropores) with two classes of macropores defined in the model.The first class incorporated tile drain-connected macropores consisting of two different modules starting from the surface or below the plough pan (32 cm depth) ending at 120 cm depth corresponding to the depth of the tile drain.The second class incorporated matrixconnected macropores consisting of three different types starting from the surface, below the plough pan (32 cm depth), and from 60 cm depth respectively, ending at the bottom of the domain (500 cm depth).It was assumed that the matrix-connected macropores deeper than 120 cm bypass the tile drain system.The macropores parameterization was based on the excavations carried out during the establishment of the monitoring programme at the field (Lindhardt et al., 2001) and the macropore studies of Nielsen et al. (2010) and Rosenbom et al. (2009).One of the objectives of testing the AM scenario was to analyse whether the preferential fast flow could be responsible for the discrepancies between measured and simulated Q and N.Leach.However, the results indicated that the AM model was unable to improve the model's overall performance (Fig. 12 and Fig. 13) and did not capture the high Q event in 2007 either (Fig. 12a).In contrast, the AM setup resulted in intensified daily events, increasing the differences between simulated and measured daily Q and hence lower KGE and NSE values (Fig. 12a compared with Fig. 2a).Furthermore, there were shortcomings in the AM model setup for the tile drain-connected macropores since the AM scenario could be considered an unrealistic modelling exercise for the 1D setup of Daisy as the distance (x) between the soil matrix and the tile drains is highly important.In the setup used for the simulations in this study, the matrix centre of the simulation was in the middle of two tile drains, which is x/2 by default in the Daisy model where x is the distance between two tile drains (here 18 m).No direct fast transport to the drains can occur over that distance.However, the defined drain-connected macropores contribute to mimicking lateral flow towards the drains through fractures in the B horizon.Hence, the fast preferential transport to drainage is not accurately described with the 1D Daisy model because it is actually a multidimensional problem.
Despite the fact that the activation of the fast flow domain (AM) did not have any effect on the total cumulative Q simulated by the model (Fig. 13a), the contribution of the AM transport significantly increased over the 10 year period (blue line, Fig. 13a).The AM contribution reached approximately 60% (blue line, Fig. 13a) in the hydrological year 2006-2007, which was 20% higher than the macropores' contribution to drainage, as reported by Larsson et al. (2007).Furthermore, it was evident that the contribution of the matrix and macropore transport to N.Leach (blue and grey lines respectively, Fig. 13b) was more or less equal from the start of the simulation in 2000 to 2004 when the high events of Q occurred.The high N.Leach event in 2004 did not bypass the drains through the deep fractures starting at 60 cm and ending at 500 cm (matrix macropores) and was simulated with AM setup as well (Fig. 12b).Although the NSE and KGE of cumulative N.leach (Fig. 13b) showed a minor improvement the accuracy of the daily resolution and events worsened (Fig. 12b; Table A.2).This indicated that the collective error decreased over the years due to the underestimation of N.Leach in 2002 and 2003.Considering the satisfactory results achieved by the initial model setup without AM (single permeability), the inclusion of macropores did not improve the simulations of Q and N.Leach for this setting of climate, agricultural management and soil.In a comparison of tile drainage discharge simulations by a single-permeability and a dual-permeability model in HYDRUS, Varvaris et al. (2018) similarly concluded that the results of the single-permeability model are satisfactory and that due to the high number of parameters needed to calibrate the dual-permeability model, the applicability of these models is questionable.Therefore, further research is essential to evaluate the macropore setup in the Daisy model when drain-connected macropores are present in the matrix, which could also benefit from an extensive parameter calibration that might result in more realistic contribution percentages by the active macropores in the setup.

Important model parameters and variable links
Performing a model sensitivity analysis for the entire parameter set is an immense task.Therefore, depending on the study and modelling objectives, specific parameters must be selected (Manevski et al., 2016).For this study, expert knowledge and the Daisy model's manual (Abrahamsen, 2020) were employed to select relevant parameters.In the current study, two main objectives were selected for the sensitivity analysis: Q and N.Leach.The analyses were conducted by methodically increasing or decreasing a single parameter within a 20% range interval maintaining the remaining parameters constant in order to identify the response of the objectives.Hydraulic conductivity of the lower boundary aquitard (k-aq) was determined to be the most sensitive parameter with Q as the objective.The k-aq determines the permeability of this layer, which directly influences the rise in groundwater level and flow to the drains.The slightest adjustments to the k-aq parameter resulted in significant nonconformities between the simulated and measured values of Q. Saturated hydraulic conductivity of the B horizon (32 to 110 cm depth) and the depth of the tile drain system were the second and third most sensitive parameters affecting the accuracy of Q predictions.Undoubtedly, P input as the main upper boundary condition had a huge effect on the prediction of Q, and the precision of the daily values were extremely important for capturing the high drainage events.While the mineral N input to the model had a direct effect on the entire N balance and N.Leach, the simulations for CrN and DM showed a high sensitivity to changes in the Fm parameter (maximum assimilation rate) in each crop during the simulation period, subsequently affecting N.Leach.
The resulted correlogram from the carried out Spearman's rankorder correlation test is presented in  simulated SWC at 25, 60 and 90 cm depths.A number of studies (De Neve and Hofman, 2002;Guntiñas et al., 2012;Dessureault-Rompré et al., 2015) have stated that volumetric or gravimetric soil water content cannot directly be used to estimate mineralisation in agricultural soil and have suggested the lack of a linear relationship between the two variables.These observations fit well with the adjustment factor F ψ m [0;1] (Fig. 3 in Hansen et al., 1991), which regulates N mineralisation at various soil water pressure potentials (expressed as pF) or equivalent levels of the soil water content, with F ψ m linearly increasing from 0.6 to 1.0 from pF 0-1.5, F ψ m equal to one at pF 1.5-2.5, and from here linearly decreasing to 0 at pF 6.5.As the soil water retention curve is strongly non-linear, a piecewise linear relationship between F ψ m and pF imposes a nonlinear relationship with SWC.

Conclusions
Despite the fact that the N fertiliser inputs during the crop rotation were according to the good management practices defined by the Danish Action Plans for the Aquatic Environment, the NO Simulated cumulative Q perfectly matched the measurements, and N.Leach was simulated satisfactorily for the 10-year period.Furthermore, SWC in 25, 60, 90 and 100 cm depth was captured acceptably by the model.The annual average water balance consisted of 650 mm of P, 464 mm of ET, 116 of deep percolation and 97 mm of Q.Moreover, the model simulated CrN and DM reasonably well, indicating realistic estimates for other N balance components such as loss through denitrification.According to the simulations, an annual average of 71% of the total N input was harvested, while losses of N included 11% to groundwater, 9% to the tile drain system, and 7% through denitrification.The main shortcomings of the model were the overestimations of N.Leach in 2004 and 2007 where the model was unable to cope with frost conditions in the soil.Furthermore, carry-over effects of the parameterised winter rape crop to the next drainage season caused the accumulation, mineralisation, and leaching of 20 and 10 kg ha −1 more NO 3 − -N than it was measured in 2004 and 2008 respectively.An assessment of the preferential flow and its effect on NO 3 − -N leaching, was carried out under a dual-permeability model setup.This assessment did not show any improvements in the Q and the daily N.Leach simulations.Furthermore, considering the great complexity of the macropore domain and the number of parameters needed for calibration, the satisfactory results achieved by the single-permeability model sufficed for explaining the measurements.The most sensitive parameter for the simulation of Q was the hydraulic conductivity of the lower boundary, whereas N. Leach showed a high sensitivity to the N input levels.
Findings of this study suggest that there is a need for improved management practices to decrease NO 3 − -N levels in the drainage water at tile-drained agricultural fields in temperate regions under conventional crop rotations, which ultimately reach the aquatic environments and coastal zones.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Design and instrumentation of the 2.33 ha Faardrup field (Lindhardt et al., 2001).P1, P2, P3 and P4 are the piezometers installed in each corner of the field.S1 and S2 include suction cups at depths of 1 m and 2 me, and TDR probes at depths of 25, 60, 90 and 110 cm.

Fig. 5 .
Fig. 5. Simulated average annual values for the inputs and outputs to the water balance (1 October-30 September) for the years 2000-2009 using the Daisy model.The input to the water balance is the precipitation (P, 650 mm), whereas the output includes evapotranspiration (ET, 464 mm), deep percolation out of the model domain (116 mm), drainage (Q, 97 mm) and storage of water in the model domain (27 mm).

Fig. 6 .
Fig. 6.NO 3 − -N concentrations at 1 and 2 m depths measured in (a) water collected from suction cups and (b) tile drainage samples in the period from 2000 to 2009.The limit for NO 3 − -N concentration in drinking water (11.3 mg L −1 ) as a red line.Grey columns show the crop rotation with winter wheat (WW), sugar beet (Sbeet), spring barley (SB), winter rape (WR) and maize (M).

Fig. 8 .
Fig. 8. (a) Measured and simulated N uptake CrN (kg N ha −1 ) and (b) measured and simulated dry matter DM (mg DM ha −1 ) in the period 2000 to 2009 after harvest of the crop in each year (RMSE = 14.12 and 1.23, respectively).

Fig. 9 .
Fig. 9. Simulated average annual values for the inputs and outputs of the N balance (1 October-30 September; kg N ha −1 ) for the years 2000-2009.The input to the N balance is 157 kg N ha −1 , changes in storage account for 26 kg N ha −1 , whereas the output includes N loss to groundwater (22 kg N ha −1 ), to the tile drain system (17 kg N ha −1 ), through denitrification (14 kg N ha −1 ) and through nitrification (4 kg N ha −1 ).
3 − -N concentration measurements was done time-proportionally until June 2004 and then changed to flowproportionally.Hence, no error could be expected in the flowcorrected NO 3 − -N leaching.
A.1b), specifically from January 2007 to March 2007 when the model was incapable of simulating the infiltration under the effect of a frost period follow by thawing.Similar to the winter rape in 2003, after the harvest of the simulated crop in July 2007, there was residual N of approximately 15 kg in the soil.Hence, as with the winter wheat in 2004, the simulations did not show any N stress for winter wheat in 2007.Therefore an explanation for the differences in N.Leach measured and simulated at the end of 2008 (Fig. A.1) could be that the model simulated abundant nitrogen in the root zone, which was transported downwards in the profile from October 2007 to March 2008 (much later compared with 2004).This resulted in a higher simulated NO 3 − -N concentration than measured at 1 m depth in the summer and autumn 2007 (Fig. A.2).The overly high simulated accumulated N in the top 100 cm of the model domain leached to the drainage system resulting in approximately 7 kg higher N.Leach simulated by the model at the end of 2008.
Fig. A.3.In agreement with the NSE and KGE metrics presented for the simulations, there was a high correlation between most of the measured and simulated water balance components (Fig. A.3).Additional correlations was observed between different variables which are described in the following paragraph.Unsurprisingly, the measured groundwater table (GWT.obs)showed high positive correlations (>0.60) with simulated daily drainage discharge (Q.sim), NO 3 − -N leaching to the drains (N.Leach.sim),soil water content at 60 and 90 cm depths (SWC60.simand SWC90.simrespectively), and NO 3 − -N leaching through deep percolation to the groundwater (N.GW.sim).Furthermore, N.Leach.obsshowed a high positive correlation with measured and simulated SWC at 90 cm depth and simulated GWT.Furthermore, simulated N mineralisation rates (Min.sim)demonstrated a positive high correlation with T and moderately high negative correlations (<−0.40) with measured and simulated GWT, Q.sim, N.Leach.sim,and measured and

Fig. 12 .Fig. 13 .
Fig. 12.(a) Measured and simulated daily Q (mm d -1 ) from 2000 to 2009 when including the macropore setup in the model, (b) measured and simulated daily N.Leach (kg N ha −1 d −1 ) from 2000 to 2009 when including the macropore setup in the model.
3 − -N concentrations in the suction cups and drainage samples from 2000 to 2009 generally exceeded the permitted level.The highest NO 3 − -N concentrations were measured from 2000 to 2002 under the winter wheat and sugar beet crop rotation and from 2005 to 2007 with maize and spring barley.

Table 1
Crop rotation, crop abbreviations in brackets, and amounts of mineral N fertiliser (NPK) applied on the Faardrup study field from 2000 to 2009.

Table 2
Physical soil properties and VG-M model parameters .

Table 3
Summary of measured and simulated annual precipitation P (mm d −1 ) and drainage Q (mm d −1 ) for hydrological years from 1 October to 30 September, and NSE (−) and KGE (−) metrics for each year.
a Year 2000 includes values from 1 January to 30 September.b Year 2009 ends on 30 September.

Table 4
Summary of measured and simulated annual cumulative N.Leach for each hydrological year (1 October to 30 September.Values of NSE and KGE are also shown.
a Year 2000 includes the values from 1 January to 30 September.b Year 2009 ends on 30 September.