Journal of Hydrology 587 (2020) 125030

Contents lists available at ScienceDirect

Journal of Hydrology
journal homepage: www.elsevier.com/locate/jhydrol

Research papers

What is the worth of drain discharge and surface runoﬀ data in hydrologicalsimulations?

T⁎

M. Turunena, , G. Gurarslanb, J. Šimůnekc, H. KoivusalodaNatural Resources Institute Finland (LUKE), Tietotie 4, 31600 Jokioinen, Finland
Department of Civil Engineering, Pamukkale University, 20160 Denizli, TurkeycDepartment of Environmental Sciences, University of California, Riverside, CA, USAd
Aalto University School of Engineering, Department of Built Environment, P.O. Box 15200, 00076 Aalto, Finlandb

A R T I C LE I N FO

A B S T R A C T

This manuscript was handled by Corrado
Corradini, Editor-in-Chief, with the assistanceof Renato Morbidelli, Associate Editor

Quantitative analyses of empirical data requirements for hydrological simulations are rare. This study aims toanalyze how a multi-objective optimization framework and information content computations aid in quantifying
ﬁeld-scale data worth in drainage studies. The results showed how a 1D numerical model and a diﬀerential
evolution algorithm performed in describing the ﬁeld water balance. The choice of the optimization target(subsurface drain discharge and surface runoﬀ) impacted the simulation results more than parameter deviations.
While the information content of surface runoﬀ data was higher than that of drain discharge, drain discharge
data contained more information on most of the soil parameters. Uncertainties related to groundwater outﬂow
data, which were not used in the optimization, were higher than those of drain discharge and surface runoﬀ. A
central weighing optimization scheme with two data types produced the best but still incomplete description of
the ﬁeld hydrology. Despite the modest model performance, the results demonstrated how the choice of empirical data and optimization strategy can lead to uncertainties in drainage simulations and how the uncertainties can be assessed. Practically, a low amount of information and a parameter sensitivity analysis canlead to a biased description of uncertainty related to such hydrological variables which are not used in the
optimization. Beneﬁts of the modeling framework were shown when assessing (1) model structure adequacywith the Pareto front analysis, (2) information content of diﬀerent data types regarding diﬀerent parameters, and(3) uncertainties related to simulating hydrological variables based on optimization against a given data type.

Keywords:
Information content
HYDRUS
Multi-objective optimizationDiﬀerential evolution
Pareto front
Identiﬁability

1. Introduction
Assessing empirical data requirements for hydrological model analyses is an important part of conducting reliable simulations, but suchrequirements are rarely quantiﬁed. Agricultural land drainage, which
has been studied for over a century (Jin and Sands, 2003), serves as aconcrete example of a hydrological problem that has been widely studied using ﬁeld measurements and simulation models. However, it is
not clear what the best way to combine measurements with simulations
is, and what the implications of the related uncertainties are. Landdrainage is currently facing new challenges due to changing climate
and rival land uses. Improved understanding of data requirements canlead to more reliable hydrological impacts assessments.
In northern conditions, drainage practices have traditionally been
implemented in agricultural ﬁelds to provide optimal traﬃcability andsoil moisture conditions for crop growth during wet spring and autumn

periods (e.g., Blann et al., 2009). However, considering environmental
impacts of agriculture (e.g., Tilman et al., 2011; Blann et al., 2009), adrainage design also oﬀers an opportunity to control environmentalload generation (Sands et al., 2015; Skaggs et al., 2005; Mendez et al.,
2004; Morrison et al., 2013; Skaggs and Chescheir, 2003). Designingsustainable land drainage is among the measures to adapt to the climatechange, which is expected to increase the frequency of mild winteroccurrence (Räisänen, 2008; Gillett et al., 2003) and environmentalloads to surface waters (e.g., Meier et al., 2012), and thus to challenge
the current drainage procedures.
While assessing hydrological processes by empirical means is ahighly laborious and costly task, combining numerical models with
empirical data provides improved means to evaluate relevant processes.
Diﬀerent numerical modeling approaches (e.g., De Schepper et al.,2017; Turunen et al., 2013; Singh et al., 2006; Boivin et al., 2006) havebeen previously developed and applied to describe hydrology of⁎Corresponding author.
E-mail addresses: mika.turunen@luke.ﬁ (M. Turunen), gurarslan@pau.edu.tr (G. Gurarslan), jsimunek@ucr.edu (J. Šimůnek),harri.koivusalo@aalto.ﬁ (H. Koivusalo).
https://doi.org/10.1016/j.jhydrol.2020.125030
Received 14 February 2020; Received in revised form 7 April 2020; Accepted 28 April 2020
Available online 04 May 20200022-1694/ © 2020 Elsevier B.V. All rights reserved.

Journal of Hydrology 587 (2020) 125030

M. Turunen, et al.

drained ﬁelds. Despite the advances in numerical modeling, empirical
data is inherently needed for model parameterization, calibration, andvalidation purposes (Silberstein, 2006).
Estimation of empirical data requirements for a robust model assessment is often not systematized but assessed a priori by expert
knowledge, which inherently involves a degree of subjectivity. The
degree of subjectivity could be decreased by exploiting methodologies
to quantify information content and worth of diﬀerent observations
with respect to the studied phenomena. For example, in earlier studies,
Rozemeijer et al. (2010) demonstrated the importance of drain discharge observations for model ﬂow route simulations. Leaf (2017)
further studied model-data relations and showed how simulations could
be applied to identify such observations which would reduce model
predictive uncertainties the most. Also, Vrugt et al. (2001) showed amethodology to quantify the information content of diﬀerent observations. The data worth of diﬀerent observations could be further analyzed by calibrating a model against diﬀerent data types one by one andassessing how accurately diﬀerent hydrological variables can be reproduced based on a chosen data type. The above quantitative methods
oﬀer approaches to identify hydrologically relevant data (e.g., Loritz
et al., 2018) and consequently to set a robust frame for the design ofappropriate ﬁeld-scale measurement campaigns for simulation studies.
Regarding drainage design studies, the data worth of diﬀerent observations has been rarely quantitatively assessed (Guillaume et al.,2019; Kikuchi, 2017). Moreover, it has been recently suggested that
information quantiﬁcation could formulate a more robust and formallyappropriate approach to provide model-based knowledge as comparedwith traditional uncertainty analyses (Nearing and Gupta, 2018;
Nearing et al., 2016). Data worth analyses can be considered to be inherently diﬀerent from widely applied sensitivity analyses (Nearing and
Gupta, 2018). The implications of data worth analyses on practice-oriented and concrete modeling cases, such as land drainage, would be
worth studying.
This study aims to analyze how numerical modeling, a multi-objective optimization framework and information content analyses (e.g.
Vrugt et al., 2001) can be applied to quantitatively describe ﬁeld-scaledata worth in the context of agricultural land drainage. The studied
data types include measured subsurface drain discharge and surface
runoﬀ time series. Drain discharge measurements are available fromnumerous research sites and they are more easily set up compared tosurface runoﬀ monitoring. This study is motivated by the need to ﬁnd
out which one of these measurements is more valuable from theviewpoint of data worth.

onwards. Thus, the hydrological data was 11 years long. The monitoredhydrological variables used in the current study were subsurface draindischarge and surface runoﬀ. Surface runoﬀ was collected with gravel
trenches (depth 0.4 m) to a subsurface drain pipe, which conveyed thewater to the measurement station. At the station, surface runoﬀ andsubsurface drain discharge were separately measured with Datawater
WS vertical helix meters (Maddalena, Povoletto, Italy) at the pipe
outlets. The outﬂow variables were measured with a time resolution of
15 min. The area of the monitored drain network was 1.7 ha, and thearea from which surface runoﬀ was collected was 1.2 ha.
Precipitation was measured on-site with a time resolution of 15 min.
Other meteorological variables (air temperature, air relative humidity,
and short-wave radiation) were available from a nearby (7 km from thesite) observatory of the Finnish Meteorological Institute. All data were
aggregated into daily time-steps for the simulations. More detailed information on the site and the data is available in Turunen et al. (2013)
and Vakkilainen et al. (2010).
2.2. Model
The ﬁeld hydrology was modeled using the HYDRUS model(Šimůnek and van Genuchten, 2008) due to its capability to describe
preferential ﬂow, generation of subsurface drainage and surface runoﬀ,
as well as other relevant hydrological processes in the snow-aﬀected
Nummela agricultural ﬁelds. HYDRUS is a modular 1D process-basedmodel, and its numerical solution is based on linear ﬁnite element
schemes. The vertical 1D model was assumed to provide a suﬃcientdescription of the ﬂat Nummela ﬁeld (e.g. Salo et al., 2015; Li et al.,2014). The time derivative is solved with a fully implicit discretization
using the Picard iterative solution scheme.
In this study, the soil water ﬂow was described with the dual-porosity module of the model, where the pore space is partitioned into
mobile (inter-aggregate) and immobile (soil matrix) regions:

∂θmo
∂ ⎡
∂h
K (h) ⎛=
+ 1⎞ ⎤ − Smo − Γ
∂t
∂x ⎢
⎝ ∂x
⎠⎥
⎣
⎦

(1)

∂θim
= −Sim + Γ
∂t

(2)
−3

where θmo and θim are the soil water content [L L ] in the mobile andimmobile regions, respectively, x is the vertical spatial coordinate [L], K(h) is the hydraulic conductivity, h is the pressure head [L], Smo and Sim
are sink terms for the two regions [T−1] and Γ is the transfer rate ofwater from the mobile to immobile regions [T−1]. Throughout the text,
the subscripts mo and im are used to refer to the mobile and immobile
pore regions, respectively. Water retention and unsaturated hydraulicconductivities were computed using the van Genuchten (1980)
schemes.
The water transfer between the two regions was described with thefollowing equation:
3

2. Materials and methods
2.1. Drained ﬁeld site and empirical data
The Nummela experimental ﬁeld is located in Southern Finland(WGS84: 60°51′59″N, 23°25′50″E). Data of the studied ﬁeld section
with an area of 1.7 ha was utilized in this study. The ﬁeld has a heavy
clay (particle size < 0.002 mm) soil with the mean clay content of≥66% (Vakkilainen et al., 2010). The ﬁeld section is ﬂat, with a slopeof < 1%. The ﬁeld was cultivated with oats and barley during themeasurement campaigns. Preferential ﬂowpaths dominate water ﬂow
in the clayey soil and allow relatively rapid drainage of the ﬁelds during
rain storm and snowmelt events (Turunen et al., 2015; Turunen et al.,2013).
The ﬁeld section (C in Turunen et al., 2013) was initially subsurface
drained in 1952 with a drain spacing of 16 m, a drain depth of 1 m, andan inner drain diameter of 0.05 m. In June 2008, new supplementarydrains were installed between the old drain lines, resulting in a drain
spacing of 8 m. Gravel was used as an envelope material, and gravel
deposits were installed every 7 m (Turunen et al., 2013; Vakkilainen
et al., 2010).
The hydrology of the ﬁeld section was monitored from May 2007

Γ = α (hmo − him)

(3)
−1 −1

where α is the ﬁrst-order mass transfer coeﬃcient [L T ]. When the
groundwater level was above the subsurface drain depth, drain discharge was described as Kroes et al. (2008):

qdr = −

hdr
γdrain(4)−1where qdr is the drain discharge [LT ], hdr is the height of the watertable [L] above the drain depth at the midpoint between the drainlines, and γ drain is the drain entrance resistance [T], which was calculated as:

γdrain =2
Ldr+ γentr
8KhBot Deq + 4KhTop hdr(5)where Ldr is the drain spacing [L], KhBot and KhTop are the horizontal
2

Journal of Hydrology 587 (2020) 125030

M. Turunen, et al.

saturated hydraulic conductivities [LT−1] below and above the drain
system, respectively, Deq is the equivalent depth [L], and γ entr is thedrain entrance resistance [T]. Deq was calculated with the method ofvan der Molen and Wesseling (1991).
The impact of groundwater ﬂow on the studied area was calculatedusing the equation of Hopmans and Stricker (1989):

qgw = −Aqh exp(−Bqh hb)

time to assess how accurately the other data type can be predictedbased on the calibration against one data type. This approach wasconsidered useful since it shows the boundaries of the possible Pareto
frontier in the model solution space (e.g. Madsen, 2003) and can thus
help to understand model-data relations and data worth (Mertens et al.,2006). In the third step, the information content of the two data types
with respect to diﬀerent parameter was quantiﬁed based on prior andposterior parameter variations (e.g. Vrugt et al., 2001; Leaf, 2017). Thiswas done to compare how accurately diﬀerent parameters can be
identiﬁed based on a data type. The combination of the above computations was assumed to give a reasonable understanding of modeldata relations and the results describe data worth in the current study(Fig. 1a). Thereafter, traditional parameter sensitivity analysis was
conducted to assess how the uncertainties related to the choice of datacompare with the parameter sensitivity. The modeling approach andthe setup are described in more detail below in Sections 2.3.2-2.3.6 andFig. 1a.

(6)
−1

−1

where qgw is the groundwater seepage ﬂux [LT ], Aqh [LT ] and Bqh[L−1] are empirical parameters, and hb [L] is the depth of the watertable below the surface of the simulated proﬁle. The seepage ﬂux is
assigned to the bottom node of the proﬁle. By simultaneously having
both drain discharge and groundwater seepage processes active, themodel version diﬀers from the standard model version of Šimůnek et al.
(2008).
The model describes snow accumulation by assuming that all precipitation is in the form of snow when the air temperature is lower than
−2 °C. Linear transition is set to control the form of snow between
temperatures of −2 and 2 °C. All precipitation is considered to be in theliquid form for temperatures higher than 2 °C. Snowmelt is describedwith a degree-day approach where the snow layer melts proportionallyto the air temperature greater than 0 °C multiplied with a degree-day
coeﬃcient D [LT−1]. Evapotranspiration from the soil inﬂuences thelayers within the root zone and is given as a precomputed time series(Turunen et al., 2015).

2.3.2. Temporal data, spatial discretization, timestepping and boundary
conditions
Since most water outﬂow at the study site occurs outside thegrowing seasons (e.g. Turunen et al., 2015), two dormant season periods were chosen for the model application. In this study, the simulations were focused on two autumn–winter-spring periods with prevailing high North Atlantic Oscillation indices, characterizing mild
winters that are presumable typical of future winter conditions underclimate warming (e.g. Gillett et al., 2003). The periods were Sep 2011 –Apr 2012 and Sep 2014 – Apr 2015 for the model calibration and validation, respectively.
Computational time steps were set to vary between 0.00001 and
0.1 days in the simulations with the dynamic time-stepping scheme of
HYDRUS (Šimůnek et al., 2008). Spatially, the ﬁeld was described using
a 2 m deep 1D proﬁle that was discretized using a 0.02 m spacing.
System-dependent boundary conditions were set at the top and bottom

2.3. Model application
2.3.1. Modeling approach
The data worth analysis was started with calibrating and validating
the model against two data types (drain discharge and surface runoﬀ) toassess the model performance for describing key hydrological processes
in the studied site (step 1 in Fig. 1a). In the second step, the model wascalibrated against one data type (drain discharge or surface runoﬀ) at aFig. 1. (a) A short description of the data worth analysis steps and (b) a schematic presentation of the structure (boundary conditions and key processes included) of
the 1D model setup.
3

Journal of Hydrology 587 (2020) 125030

M. Turunen, et al.

Vi, G + 1 = PR1, G + F (PR2, G − PR3, G ), R1 ≠ R2 ≠ R3 ≠ iof the proﬁle (Šimůnek et al., 2008). At the top, the atmospheric condition was controlled by precipitation, evaporation and h at the soil
surface. At the bottom, the boundary was described by the sum of Eqs.
(4) and (6) (drain discharge and groundwater seepage, respectively). A
schematic presentation of the boundary conditions and key processes
included in the model is shown in Fig. 1b.

where i, R1, R2, and R3 denote randomly selected vectors among the
population, F is the mutation factor and G is the generation number.
The parent vector Pi,G and the mutant vector Vi,G+1 are merged into a
trial vector Ui,G+1 parameter-wise as follows:

Vi, k, G + 1, if R ≤ CR or k = Rk
Ui, k, G + 1 = ⎧otherwise⎨⎩ Pi, k, G ,2.3.3. Multi-objective optimization
The model was calibrated against multiple data types using thediﬀerential evolution algorithm. The data type term refers to measurements of diﬀerent hydrological phenomena. The two data types inthis study included cumulative surface runoﬀ and drain discharge.
Regarding the objective function, the widely used root mean squareerror (RMSE) was chosen to describe the diﬀerence between the simulated and observed time series (fj [L]). Individual objective functions
were formed by compensating diﬀerences in the magnitude of the twoseries by normalizing the RMSE values by the standard deviation of thedata (nRMSE) and further by transforming the distance of the nRMSEs
from the origin near the optimum (e.g., Köhne et al., 2011; Mertens
et al., 2006; Madsen 2003):

Fj =

fj
σj

+ φj , j = 1, 2

⎜

(11)

where Uik,G+1 is the resulting trial vector parameter, k is the parameter
index in a vector, R ∈ [0,1] is a random number, CR ∈ [0,1] is the
crossover parameter, and Rk is a randomly selected index. A new R is
generated for each k, but Rk is common for each k of the trial vector.
Thus, Rk ensures that at least one parameter is obtained from the mutant vector. The trial vector Ui,G+1 replaces the parent vector Pi, G if theobjective function value of the trial vector equals or is lower than thevalue of the parent. In the current study, the algorithm was set to stop
when the standard deviation of the population objective function values
was lower than 1% of the mean population objective function value.
The parameters F and CR were set to the value of 0.8 following
Gurarslan and Karahan (2015). Diﬀerent population numbers (20, 40,
and 60) were tested to address their impact on the optimization results.
The population number had practically no impact on the results, andthe population number of 40 was used in the optimizations. Initialization of the population was performed using the Latin Hypercube sampling. Non-convergent simulations were given a very high objectivefunction value. Also, since running a successful simulation typicallytook only a couple of seconds, simulations taking more than 30 s wereaborted and considered non-convergent, following Brunetti et al.
(2016).

(7)

where Fj is the individual objective function, σj is the standard deviationof the data [L] associated with f j , φj is the transformation constant, andj denotes diﬀerent objectives (e.g., drain ﬂow and surface runoﬀ). φj
was calculated as follows (e.g., Köhne et al., 2011; Mertens et al.,2006):fjf f
φj = max ⎛ 1 , 2 ⎞ − , j = 1, 2σσσj
⎝ 1 2⎠

(10)

⎟

(8)

2.3.5. Model setup and calibrated parameters
The spatial domain of the model was set to consist of ﬁve soil
horizons with distinct soil properties: 0–0.2, 0.2–0.4, 0.4–0.6, 0.6–1.0,
and 1.0–2.0 m. The soil hydraulic parameters for the soil matrix (immobile zone) were derived from Turunen et al. (2015), as listed inTable 1.
Based on previous studies in clayey ﬁelds in southern Finland, it was
assumed that in addition to the soil hydraulic parameters in Table 1,
water ﬂow in the ﬁeld was mostly controlled by the hydraulic properties of the macropore system, the parameters of the subsurface drainage
system, and the parameters controlling groundwater outﬂow from theﬁeld (Warsta et al., 2013; Turunen et al., 2013). Eight parameters werechosen to describe these features, and plausible parameter value ranges
were identiﬁed from the literature. The saturated hydraulic conductivity (Ksat) of the mobile zone was computed as a product of aconductivity multiplier (Ksat,mul [L T−1]) and θs_mo (Table 1) following
Warsta et al. (2013) and Larsbo and Jarvis (2003). The calibration
parameters and their ranges are given in Table 2. The van Genuchten
shape parameters for the mobile zone were adopted from Gärdenäs
et al. (2006).

The transformation constant keeps each individual objective function the same distance from the origin and thus balances their contributions in the aggregated objective function. The actual aggregatedobjective function (Fagg ) was then formed by aggregating the individual
objective functions following the central weighing scheme (Ocw)(Mertens et al., 2004; Madsen 2003):

Fagg = wF1 + (1 − w ) F2

(9)where w is the weighting coeﬃcient, which is given a value of 0.5 in thecentral weighing scheme.
2.3.4. Diﬀerential evolution
The optimization was conducted with the metaheuristic diﬀerential
evolution algorithm (DEA), which is considered to be an eﬀective
global optimization method (e.g., Storn and Price, 1997) and has previously been successfully implemented in several hydrological modeling studies (e.g., Gurarslan and Karahan, 2015; Ayvaz and Gurarslan,
2019; Karterakis et al., 2007).
The DEA starts by setting a number of candidate solutions in amodel parameter space. The algorithm considers each population
member as a vector of parameters and evaluates the ﬁtness of the population members with a given objective function (Eq. (9)). The search
converges towards the global minimum by stochastically combining theparameters of population members by adding weighed diﬀerences between two population vectors to a third vector. If the resulting vector
has a lower objective function value than a chosen older (parent)
vector, the new vector replaces the old vector. Thus, the population
members self-organize towards the global minimum. There are numerous diﬀerent approaches to apply the DEA (e.g., Storn and Price,1997). In the current study, the so-called DE/rand/1/bin (Storn andPrice, 1997) strategy is applied, following Gurarslan and Karahan(2015). The approach is presented shortly below.
In the algorithm, a mutant vector Vi,G+1 is created for each population parent vector Pi,G during each generation:

Table 1The soil hydraulic parameters for diﬀerent soil depths.
Parameter
Depth [m]θs_im [–]

θr_im [–]

αim [m−1]

nim [–]

θs_mo [–]

0.0–0.2 m
0.2–0.4 m
0.4–0.6 m
0.6–1.0 m
1.0–2.0 m

0.529
0.539
0.534
0.534
0.534

0.1
0.1
0.1
0.1
0.1

2.886
0.412
0.679
0.679
0.679

1.149
1.161
1.127
1.127
1.127

0.03360.00230.0034

θs_im, θr_im, αim, and nim are parameters of the van Genuchten (1980) model forthe immobile pore system (θs_im and θr_im are saturated and residual watercontents, and αim and nim are the shape parameters).
θs_mo is a parameter for the mobile pore system.
4

Journal of Hydrology 587 (2020) 125030

M. Turunen, et al.

Table 2
Calibrated parameters and ranges used in model optimization.
Parameter

SymbolKsat multiplierθs_mo (depth 0.6–1.0 m)θs_mo (depth 1.0–2.0 m)
Tortuosity
Mass transfer coeﬀ. in Eq. (3)Parameter in Eq. (6)Parameter in Eq. (6)Drainage resistance

Range

Reference value source−1

Ksat,mulθs_mo_1
θs_mo_2lαAqh
Bqh
γ entr120…1920 [m d ]0.0001…0.0031 [m3 m−3]0.0001…0.0031 [m3 m−3]
−1.0…1.0 [–]0.001…0.1 [m−1d−1]0.00001 …0.1 [md−1]0.001 …10.0 [m−1]10…250 [d]Salo et al. (2017)
Turunen et al. (2013)
Turunen et al. (2013)
Schaap and Leij (2000)
Gärdenäs et al. (2006)
Hopmans and Stricker (1989)
Hopmans and Stricker (1989)
Kroes et al. (2008)

2.3.6. Data worth and sensitivity
To assess uncertainties related to the possibly Pareto-optimized solution, the model was also optimized separately against the discharge
data type (Odr) and the runoﬀ data type (Or), i.e., by setting φ to 0 (Eq.
(7)) and w to 0 or 1 (Eq. (9)) when optimizing against drainage orsurface runoﬀ data, respectively. The purpose of the analysis was toexplore the Pareto front extremes and, consequently, to analyze how
weighing of individual objective functions aﬀected the solution (e.g.,
Madsen 2003).
To further analyze the role of two diﬀerent calibration data types in
the model optimization and model-data interactions, their information
content was analyzed similarly as by Vrugt et al. (2001):

ICm (p) = 1 −σ (p)mσ (p)(12)

where ICm(p) is information content of measurement m regarding
parameter p, σ(p)m is the standard deviation of parameter p in the accepted (behavioral) simulations with respect to measurement m andσ(p) is the standard deviation of the parameter in the initial parameter
ranges. Thus, an ICm(p) value close to one indicates that the information content of measurement m regarding parameter p is high. An ICm(p) value near zero indicates a low information content. Information
contents were analyzed with a simulation test where the parameter
ranges were set as in Table 2. Then, 80,000 simulations were conductedusing the Latin Hypercube approach and such simulations, which had
the nRMSE < 1.0, were considered behavioral in the information
content calculation (Eq. (12)).
The sensitivity of the solution optimized with Ocw (the central
weighing scheme) to parameter deviations was assessed with a global
sensitivity analysis where the calibrated parameters were alteredby ± 10%. The model realizations were conducted using the Latin
Hypercube Monte Carlo sampling. Following Brunetti et al. (2016),
5000 realizations per parameter were assumed to be suﬃcient to reasonably explore the model parameter space. Thus, the analysis was
conducted with 40,000 (8 parameters and 5000 realizations) modelruns.

Fig. 2. The dominated and optimized solutions in the model solution space.
Table 3The optimized parameter values with the diﬀerent optimization objectives.
Parameter

Central weighing
schemeDrain discharge
objective

Surface runoﬀ
objective

Unit

Ksat,mulθs_mo_1
θs_mo_2lαAqh
Bqh
γ entr125.9
0.00220.00020.53
0.053
0.0002
0.37
244.6

122.0
0.00290.00050.71
0.039
0.0026
1.75
78.0

132.1
0.00080.00020.98
0.075
0.00128.76
239.8[m d−1][m3 m−3][m3 m−3][–][m−1d−1][m d−1][m−1][d]outﬂow (37 mm) and soil water storage change (19 mm). Thus, draindischarge was clearly the dominant water balance component. During
the simulation, the groundwater table depth ﬂuctuated between the soil
surface and the depth of 1.4 m. The validation results were qualitatively
similar to the calibration results, but the model overestimated draindischarge and underestimated surface runoﬀ (Fig. 3).

3. Results
3.1. Calibration and validation
The DEA optimization resulted in nRMSE of 0.19 and 0.18 for draindischarge and surface runoﬀ, respectively, after 81,831 model realizations (Fig. 2). The resulting model parameters are given in Table 3. The
optimum is located at the sharp edge of the model solution space (bluecircle in Fig. 2). The Pareto front was found to be unsymmetrical, withconsiderably longer Euclidean distance between the optimum of Ocw
and Odr (dark green circle in Fig. 2) than between Ocw and Or (lightgreen circle in Fig. 2).
The measured and simulated cumulative drain discharge and surface runoﬀ for the calibration and validation periods are illustrated inFig. 3 with the corresponding evaluation criteria. Other key waterbalance components were evapotranspiration (71 mm), groundwater

3.2. Sensitivity and information content analyses
Simulation results can be considered to be rather conservative withrespect to parameter deviations, since in qualitative terms, they do notcause pronounced changes in the distribution of the two water balancecomponents (Fig. 4). Simulation results were found to be less sensitive
to parameter deviations than to the choice of the optimization objective(Fig. 4), i.e., to the choice between the extremes of the Pareto front, asshown in Fig. 2. Regarding the impact of the optimization approach(i.e., Odr or Ocw) on outﬂow components, the accumulated surface
runoﬀ for Odr diﬀered by −63 mm from the one for Ocw. Cumulativedrain discharge for Or diﬀered by 28 mm from the one for Ocw.
5

Journal of Hydrology 587 (2020) 125030

M. Turunen, et al.

4. Discussion
Model optimization with the diﬀerential evolution algorithm, central weighing scheme and the two objectives resulted in a reasonablereproduction of the measured variables, but there were periods with aclear mismatch, more often during validation than calibration. The
mismatch partly indicates epistemic uncertainties, such as meteorological data inconsistent with hydrological measurements or limited
knowledge of system dynamics. While epistemic uncertainties are difﬁcult to quantify (e.g., Beven 2016), data worth of the partly incompletely measured variables is the remaining question which ouranalysis aimed to answer.
The sharp edge of the Pareto front (Fig. 2) indicated that optimizations with two objectives did not induce a signiﬁcant decrease in the
goodness of a single-objective optimum. Such an ability to simultaneously optimize two objectives is an indication of an adequate model
structure regarding simulations of studied phenomena (Mertens et al.,2006). The Euclidean distance between the optimum and the origin
indicates the magnitude of modeling and measurement error, which canbe considered small compared with the variance of the data (Fig. 2).
The applied multi-objective framework can be used to partly demonstrate the adequacy of a model structure. The value of the approach ishighlighted by the fact that such an assessment is practically impossible
with a traditional model calibration–validation framework.
Optimization with only one objective led to a reasonable descriptionof that variable, but to uncertainties regarding other studied variables.
Rozemeijer et al. (2010) showed how drain discharge data could be
more important than groundwater level data when simulating wateroutﬂow pathway contributions. Our results further pointed out thatwhen optimizing surface runoﬀ (Or) it was possible to estimate draindischarge slightly more accurately than it was to estimate surface runoﬀ
when optimizing drain discharge (Odr) (Fig. 4). On the other hand,
optimization with the surface runoﬀ objective (Or) resulted in an unreasonably low (8 mm) estimate of groundwater outﬂow, as comparedwith earlier studies on the site (Turunen et al., 2015; Turunen et al.,2013). This distinction was also reﬂected in the overestimation of simulated drain discharge in the Or scheme in the current study. Only thecentral weighing scheme resulted in a reasonable description of allconsidered outﬂow ﬂuxes. Regarding model outputs, relatively highuncertainties in data types not used in model calibration were previously reported by Højberg and Refsgaard (2005) in a groundwater
simulation study. Intensive ﬁeld-scale data sets with various hydrological data types could be analyzed with the current approach to show
how accurately diﬀerent variables can be predicted based on anothervariable. In other words, the approach can practically help quantify theworth of diﬀerent data and, consequently, to identify which data is

Fig. 3. Measured and simulated drain discharge and surface runoﬀ during (a)calibration and (b) validation.

Similarly, diﬀerences in simulated groundwater outﬂow were considerable for diﬀerent optimization approaches, i.e., groundwater outﬂow was 37 mm for Ocw, 126 mm for Odr, and 8 mm for Or. Diﬀerencesin the parameters controlling the amount of groundwater outﬂow (i.e.,
Aqh and Bqh in Table 3) explained diﬀerences in the amount of simulated
groundwater outﬂow.
Regarding the information content of the diﬀerent observations, themean ICm(p) values for all parameters were 0.05 and 0.19 for Odr andOr, respectively. The number of converged realizations with the nRMSE
below the threshold of 1.0 was 43,603 and 540 for Odr and Or, respectively. The Monte Carlo process was considered to converge basedon the mean information content and outﬂow values shown in Fig. 5.
Regarding individual parameters, the ICm(p) values are shown in
Table 4. The results showed that the surface runoﬀ data type contained
more information for θs_mo_1, Aqh, and γ entr than drain discharge data
type, while the drain discharge data contained more information forremaining parameters than surface runoﬀ data. The information content was the highest regarding Aqh and γ entr when the surface runoﬀ
objective was used, and regarding Bqh when the drain discharge objective was used. Information contents were lowest regarding Ksat,mul,θs_mo_2, I, and α. Note that parameter ranges of accepted simulations
were relatively wide (Table 4 vs. Table 2).

Fig. 4. Simulation results for cumulative ﬂuxes (drain discharge and surface runoﬀ) and their deviations due to the choice of the model optimization strategy and thesensitivity analysis.
6Journal of Hydrology 587 (2020) 125030

M. Turunen, et al.

Fig. 5. The convergence of mean cumulative outﬂow and information content for (a-b) drain discharge and (c-d) surface runoﬀ during the Monte Carlo process. Redlines denote the mean value and blue lines the value of an individual realization. (For interpretation of the references to colour in this ﬁgure legend, the reader isreferred to the web version of this article.)

the soil entered the trench in addition to the runoﬀ from the ﬁeld
surface. Thus, the information content analysis may include such
epistemic uncertainties.
The computational concept could be further applied to study information contents during a range of diﬀerent periods, as the nature ofdata worth can be partially dynamic (Loritz et al., 2018; Vrugt et al.,2001). For example, during growing seasons when groundwater outﬂow dominates the total outﬂow (e.g., Turunen et al., 2013) it is likelythat data types describing subsurface hydrological observations havehigher information content than for example surface runoﬀ. This points
out the importance of considering the timing of monitoring periods
when conducting a data worth analysis.
The information content analysis also showed how much diﬀerent
data types contained information about diﬀerent parameters. The
highest diﬀerence between the two data types was found in the parameters controlling groundwater outﬂow and drainage resistance.
Unexpectedly, the surface runoﬀ data constrained the range of thedrainage resistance values considerably more than the drain discharge

relevant for a hydrological analysis.
The simulations further showed how a model framework could be
set up to quantify diﬀerences in information contents of diﬀerent datatypes. The information content of surface runoﬀ observations washigher than that of drain discharge observations. However, regarding
parameters describing soil hydraulic properties (Ksat,mul, θs_mo_1, θs_mo_2,
l, and α) rather than sink term properties, the information content ofdrain discharge was typically higher than that of surface runoﬀ(Table 4). It is noteworthy that the absolute values of information
content were aﬀected by the choice of parameter ranges, which also
typically involve subjectivity despite their physical meaning (e.g.,
Beven and Binley 1992). However, the information content analysis
facilitated a quantitative comparison between diﬀerent data types.
Practically, the approach can help a modeler to understand how much
uncertainty is reﬂected in the model parameters when the model is
optimized against a chosen data type. It is also possible that surface
runoﬀ simulations and surface runoﬀ observations (measured with a
0.4 m deep trench) were not entirely commensurable if seepage from

Table 4The computed information contents of the diﬀerent parameters and the corresponding parameter value ranges.
Information content

Parameter value range

Parameter

Drain discharge objective

Surface runoﬀ objective

Drain discharge objective

Surface runoﬀ objective

Unit

Ksat,mulθs_mo_1
θs_mo_2lαAqh
Bqh
γ entr

0.01
0.01
0.02
0.01
0.01
< 0.000.28
0.04< 0.000.05< 0.00< 0.00< 0.00
0.94
0.04
0.66

120.1…1920.0
0.0001…0.0031
0.0001…0.0031−1.0…1.0
0.001…0.1
0.00001…0.1
0.003…10.0
10.0…250.0

120.1…1919.7
0.0001…0.0031
0.0001…0.0031−1.0…1.0
0.001…0.1
0.00002…0.02
0.003…10.0
134.9…250.0[m d−1][m3 m−3][m3 m−3][–][m−1d−1][md−1]
[m−1]
[d]

7

Journal of Hydrology 587 (2020) 125030

M. Turunen, et al.

types resulted in the best but still incomplete description of the ﬁeld
hydrology, pointing out epistemic uncertainties. In a wider context, thequantitative modeling framework could be applied to map the data
worth of a wide range of diﬀerent data types to evaluate what types ofdata are required to conduct reliable assessments of chosen hydrological phenomena.

data, since the saturation-excess mechanism triggering surface runoﬀ
generation requires high subsurface moisture conditions, which arelargely controlled by drainage. This mechanism also pointed out the
coupling between the surface and subsurface domains and partly exempliﬁed why both studied data types contain information about eachother.
The parameter ranges of the single-objective optimizations (Table 4)
were rather wide, which implied considerable non-uniqueness in thosesolutions. The non-uniqueness and possible noise in empirical dataexplain why it can be diﬃcult to identify a plausible model with only
one objective. Due to possible small inaccuracies in input or calibration
data it can be practically impossible to identify the most plausiblemodel from many rather similarly performing model realizations(Fig. 2). However, increasing the number of data types in the simulation
can narrow down the model solution space (Fig. 2). Model ensemble
approaches can be argued to be often a necessity for reliable model
assessments for decision making (e.g. Ferré, 2017). Information quantiﬁcation of a data type can help identify the limits of model predictions, which can be considered to be an inherently diﬀerent approachthan a parameter sensitivity analysis (e.g., Nearing and Gupta, 2018).
Our results provided a practical demonstration of how relatively
low amounts of information and a parameter sensitivity analysis canlead to a biased description of uncertainty related to such hydrologicalvariables which are not used in the optimization (Fig. 4). Informationon the studied parameters could be more eﬃciently extracted from
temporal data series by identifying such data subsets containing themost information for one of the parameters (Vrugt et al., 2001). While
this study assessed model uncertainty on the ﬁeld-scale, it is suggested
that studying the data worth of widely available hydrological satellite
data (e.g., Xu et al., 2014) combined with simple terrestrial observations could facilitate the assessment of the reliability of drainage-related modelling on a larger scale.
Overall, our data worth analysis showed how the modeling framework could be used to assess (1) adequacy of the model structure with
the Pareto front analysis, (2) information content of diﬀerent data types
regarding diﬀerent parameters, and (3) uncertainties related to simulating hydrological variables based on optimization against a given data
type. According to Kikuchi (2017), such a quantitative approach can beconsidered as a formal data worth analysis as compared with more
widely conducted informal analyses, which are based on expert judgement. As pointed out by Kikuchi (2017) and Guillaume et al. (2019),
data worth analyses are relatively rarely conducted despite their potential beneﬁts. Our study demonstrated the beneﬁts in the context ofdrainage and ﬁeld-scale studies.

CRediT authorship contribution statement
M. Turunen: Conceptualization, Methodology, Software, Formal
analysis, Visualization, Validation, Writing - original draft. G.
Gurarslan: Methodology, Formal analysis, Writing - review & editing.
J. Šimůnek: Methodology, Software, Writing - review & editing. H.
Koivusalo: Conceptualization, Funding acquisition, Writing - review &editing.
Declaration of Competing Interest
The authors declare that they have no known competing ﬁnancial
interests or personal relationships that could have appeared to inﬂuence the work reported in this paper.
Acknowledgements
This study was funded by Maa- ja vesitekniikan tuki ry and Natural
Resources Institute Finland. Finnish Field Drainage Association is acknowledged for sharing the hydrological and soil data. Merja Myllys is
acknowledged for support with the empirical data. The authors also
wish to acknowledge CSC – IT Center for Science, Finland, for computational resources.
References
Ayvaz, M.T., Gurarslan, G., 2019. A hybrid optimization approach for parameter estimation of conﬁned and leaky conﬁned aquifers. Water Supply 19, 2359–2366.
Beven, K., Binley, A., 1992. The future of distributed models: model calibration anduncertainty prediction. Hydrol. Process. 6, 279–298.
Blann, K.L., Anderson, J.L., Sands, G.R., Vondracek, B., 2009. Eﬀects of agriculturaldrainage on aquatic ecosystems: a review. Critical Reviews in Environmental Scienceand Technology 39, 909–1001.
Boivin, A., Šimůnek, J., Schiavon, M., van Genuchten, M.T., 2006. Comparison of pesticide transport processes in three tile-drained ﬁeld soils using HYDRUS-2D. Vadose
Zone J. 5, 838–849.
Brunetti, G., Šimůnek, J., Piro, P., 2016. A comprehensive numerical analysis of the hydraulic behavior of a permeable pavement. J. Hydrol. 540, 1146–1161.
De Schepper, G., Therrien, R., Refsgaard, J.C., He, X., Kjaergaard, C., Iversen, B.V., 2017.
Simulating seasonal variations of tile drainage discharge in an agricultural catchment. Water Resour. Res. 53, 3896–3920.
Beven, K., 2016. Facets of uncertainty: epistemic uncertainty, non-stationarity, likelihood, hypothesis testing, and communication. Hydrol. Sci. J. 61, 1652–1665.
Ferré, T.P., 2017. Revisiting the relationship between data, models, and decision-making.
Groundwater 55, 604–614.
Gillett, N.P., Graf, H.F., and Osborn, T.J., 2003. Climate change and the North Atlantic
oscillation. Geophysical monograph – American Geophysical Union 134, 193-210.
Guillaume, J.H., Jakeman, J.D., Marsili-Libelli, S., Asher, M., Brunner, P., Croke, B., Hill,
M.C., Jakeman, A.J., Keesman, K.J., Razavi, S., Stigter, J.D., 2019. Introductory
overview of identiﬁability analysis: A guide to evaluating whether you have the righttype of data for your modeling purpose. Environ. Modell. Software 119, 418–432.
Gurarslan, G., Karahan, H., 2015. Solving inverse problems of groundwater-pollutionsource identiﬁcation using a diﬀerential evolution algorithm. Hydrogeol. J. 23,
1109–1119.
Gärdenäs, A.I., Šimůnek, J., Jarvis, N., van Genuchten, M.T., 2006. Two-dimensionalmodelling of preferential water ﬂow and pesticide transport from a tile-drained ﬁeld.
J. Hydrol. 329, 647–660.
Højberg, A.L., Refsgaard, J.C., 2005. Model uncertainty–parameter uncertainty versus
conceptual models. Water Sci. Technol. 52, 177–186.
Hopmans, J.W., Stricker, J.N.M., 1989. Stochastic analysis of soil water regime in awatershed. J. Hydrol. 105, 57–84.
Jin, C.X., Sands, G.R., 2003. The long-term ﬁeld-scale hydrology of subsurface drainage
systems in a cold climate. Trans. ASAE 46, 1011–1021.
Karterakis, S.M., Karatzas, G.P., Nikolos, I.K., Papadopoulou, M.P., 2007. Application of
linear programming and diﬀerential evolutionary optimization methodologies for thesolution of coastal subsurface water management problems subject to environmental
criteria. J. Hydrol. 342, 270–282.
Kikuchi, C., 2017. Toward increased use of data worth analyses in groundwater studies.

5. Conclusions
The numerical simulations combined with the ﬁeld-scale hydrological data demonstrated how the multi-objective optimization framework, information content computations and the Pareto front analysis could be beneﬁcial for drainage studies and designs of empiricalmonitoring campaigns. Their beneﬁts were demonstrated when assessing (1) model structure adequacy with the Pareto front analysis, (2)information content of diﬀerent data types regarding diﬀerent parameters, and (3) uncertainties related to simulating hydrological variables based on optimization against a given data type. It was shownhow the choice of the empirical data and optimization strategy had alarger impact on simulation results than parameter deviations in themodel sensitivity analysis. The information content of the surface
runoﬀ data was higher than that of subsurface drain discharge. Whendiﬀerent optimization strategies were compared, the highest uncertainties were related to those data, which were not used in any of the
optimization approaches. Practically, a low amount of information anda parameter sensitivity analysis can lead to a biased description of
uncertainty related to such hydrological variables which are not used in
the optimization. The central weighing scheme with two available data
8

Journal of Hydrology 587 (2020) 125030

M. Turunen, et al.

Salo, H., Warsta, L., Turunen, M., Paasonen-Kivekäs, M., Nurminen, J., Koivusalo, H.,
2015. Development and application of a solute transport model to describe ﬁeld-scale
nitrogen processes during autumn rains. Acta Agricul. Scand., Sect. B—Soil & Plant
Sci. 65 (sup1), 30–43.
Salo, H., Warsta, L., Turunen, M., Nurminen, J., Myllys, M., Paasonen-Kivekäs, M.,
Alakukku, L., Koivusalo, H., 2017. Simulating 3-D water ﬂow in subsurface drain
trenches and surrounding soils in a clayey ﬁeld. Soil Tillage Res. 168, 20–32.
Silberstein, R.P., 2006. Hydrological models are so good, do we still need data? Environ.
Modell. Software 21, 1340–1352.
Šimůnek, J., M. Šejna, H. Saito, M. Sakai, and M. Th. van Genuchten, The HYDRUS-1D
Software Package for Simulating the Movement of Water, Heat, and Multiple Solutes
in Variably Saturated Media, Version 4.08, HYDRUS Software Series 3, Department ofEnvironmental Sciences, University of California Riverside, Riverside, California,
USA, pp. 330, 2008.
Šimůnek, J., van Genuchten, M.T., 2008. Modeling nonequilibrium ﬂow and transport
processes using HYDRUS. Vadose Zone J. 7, 782–797.
Singh, R., Helmers, M.J., Qi, Z., 2006. Calibration and validation of DRAINMOD to designsubsurface drainage systems for Iowa's tile landscapes. Agric. Water Manag. 85,
221–232.
Skaggs, R.W., Chescheir III, G.M., 2003. Eﬀects of subsurface drain depth on nitrogen
losses from drained lands. Trans. ASAE 46, 237.
Skaggs, R.W., Youssef, M.A., Chescheir, G.M., Gilliam, J.W., 2005. Eﬀect of drainage
intensity on nitrogen losses from drained lands. Trans. ASAE 48, 2169–2177.
Storn, R., Price, K., 1997. Diﬀerential evolution: a simple and eﬃcient heuristic strategy
for global optimization over continuous spaces. J. Global Optimimiz. 11, 341–359.
Tilman, D., Balzer, C., Hill, J., Befort, B.L., 2011. Global food demand and the sustainable
intensiﬁcation of agriculture. Proc. Natl. Acad. Sci. 108, 20260–20264.
Turunen, M., Warsta, L., Paasonen-Kivekäs, M., Nurminen, J., Myllys, M., Alakukku, L.,Äijö, H., Puustinen, M., Koivusalo, H., 2013. Modeling water balance and eﬀects ofdiﬀerent subsurface drainage methods on water outﬂow components in a clayey
agricultural ﬁeld in boreal conditions. Agric. Water Manag. 121, 135–148.
Turunen, M., Warsta, L., Paasonen-Kivekäs, M., Nurminen, J., Koivusalo, H., 2015.
Simulating water balance and evapotranspiration in a subsurface drained clayey
agricultural ﬁeld in high-latitude conditions. Acta Agriculturae Scandinavica, Section
B—Soil & Plant. Science 65 (sup1), 44–57.
Vakkilainen, P., Alakukku, L., Koskiaho, J., Myllys, M., Nurminen, J., Paasonen-Kivekäs,
M., Peltomaa, R., Puustinen, M., and Äijö, H., 2010. Pellon vesitalouden optimointi –Loppuraportti 2010. Salaojituksen tutkimusyhdistys ry:n tiedote 30. Finnish Field
Drainage Association, Helsinki, Finland, 114 pp.
van der Molen, W.H., Wesseling, J., 1991. A solution in closed form and a series solution
to replace the tables for the thickness of the equivalent layer in Hooghoudt's drain
spacing equation. Agric. Water Manag. 19, 1–16.
Vrugt, J.A., Bouten, W., Weerts, A.H., 2001. Information content of data for identifying
soil hydraulic parameters from outﬂow experiments. Soil Sci. Soc. Am. J. 65, 19–27.
Warsta, L., Taskinen, A., Koivusalo, H., Paasonen-Kivekäs, M., Karvonen, T., 2013.
Modelling soil erosion in a clayey, subsurface-drained agricultural ﬁeld with a threedimensional FLUSH model. J. Hydrol. 498, 132–143.
Xu, X., Li, J., Tolson, B.A., 2014. Progress in integrating remote sensing data and hydrologic modeling. Prog. Phys. Geogr. 38, 464–498.

Groundwater 55, 670–673.
Köhne, J.M., Wöhling, T., Pot, V., Benoit, P., Leguédois, S., Le Bissonnais, Y., Šimůnek, J.,
2011. Coupled simulation of surface runoﬀ and soil water ﬂow using multi-objective
parameter estimation. J. Hydrol. 403, 141–156.
Kroes, J.G., Van Dam, J.C., Groenendijk, P., Hendriks, R.F.A., and Jacobs, C.M.J., 2008.
SWAP Version 3.2 – Theory 753 Description and User Manual. Alterra-report 1649.
Wageningen, The Netherlands, 262 pp.
Larsbo, M., and Jarvis, N., 2003. MACRO 5.0. A model of water ﬂow and solute transport
in macroporous soil. Technical description, Swedish University of Agricultural
Sciences, 40 p.
Leaf, A.T., 2017. Using models to identify the best data: An example from northern
Wisconsin. Groundwater 55, 641–645.
Li, Y., Šimůnek, J., Jing, L., Zhang, Z., Ni, L., 2014. Evaluation of water movement andwater losses in a direct-seeded-rice ﬁeld experiment using Hydrus-1D. Agric. Water
Manag. 142, 38–46.
Loritz, R., Gupta, H., Jackisch, C., Westhoﬀ, M., Kleidon, A., Ehret, U., Zehe, E., 2018. On
the dynamic nature of hydrological similarity. Hydrol. Earth Syst. Sci. 22,
3663–3684.
Madsen, H., 2003. Parameter estimation in distributed hydrological catchment modellingusing automatic calibration with multiple objectives. Adv. Water Resour. 26,
205–216.
Meier, H.E.M., Hordoir, R., Andersson, H.C., Dieterich, C., Eilola, K., Gustafsson, B.G.,
Höglund, A., Schimanke, S., 2012. Modeling the combined impact of changing climate and changing nutrient loads on the Baltic Sea environment in an ensemble oftransient simulations for 1961–2099. Clim. Dyn. 39, 2421–2441.
Mendez, A., Sands, G.R., Basin, B., Jin, C.X., Wotzka, P.J., 2004. Simulating the impact ofdrainage design in a cold climate with ADAPT1. JAWRA J. Am. Water Resour. Assoc.
40, 385–400.
Mertens, J., Madsen, H., Feyen, L., Jacques, D., Feyen, J., 2004. Including prior information in the estimation of eﬀective soil parameters in unsaturated zone modelling. J. Hydrol. 294, 251–269.
Mertens, J., Stenger, R., Barkle, G.F., 2006. Multiobjective inverse modeling for soil
parameter estimation and model veriﬁcation. Vadose Zone J. 5, 917–933.
Morrison, J., Madramootoo, C.A., Chikhaoui, M., 2013. Modeling the inﬂuence of tile
drainage ﬂow and tile spacing on phosphorus losses from two agricultural ﬁelds insouthern Québec. Water Quality Res. J. 48, 279–293.
Nearing, G.S., Tian, Y., Gupta, H.V., Clark, M.P., Harrison, K.W., Weijs, S.V., 2016. A
philosophical basis for hydrological uncertainty. Hydrol. Sci. J. 61, 1666–1678.
Nearing, G.S., Gupta, H.V., 2018. Ensembles vs. information theory: supporting science
under uncertainty. Front. Earth Sci. 12, 653–660.
Räisänen, J., 2008. Warmer climate: less or more snow? Clim. Dyn. 30, 307–319.
Rozemeijer, J.C., van der Velde, Y., McLaren, R.G., van Geer, F.C., Broers, H.P., Bierkens,
M.F.P., 2010. Integrated modeling of groundwater-surface water interactions in atile-drained agricultural ﬁeld: The importance of directly measured ﬂow route contributions. Water Resour. Res. 46.
Schaap, M.G., Leij, F.J., 2000. Improved prediction of unsaturated hydraulic conductivity
with the Mualem-van Genuchten model. Soil Sci. Soc. Am. J. 64, 843–851.
Sands, G.R., Canelon, D., Talbot, M., 2015. Developing optimum subsurface drainage
design procedures. Acta Agriculturae Scandinavica, Section B—Soil & Plant. Science
65 (sup1), 121–127.

9