Article

Investigating Atrazine Concentrations in the
Zwischenscholle Aquifer Using MODFLOW with the
HYDRUS-1D Package and MT3DMS
Sahila Beegum 1,*, Jan Vanderborght 2, Jiří Šimůnek 3, Michael Herbst 2, K. P. Sudheer 1,4
and Indumathi M Nambi 1
Department of Civil Engineering, Indian Institute of Technology Madras, Chennai 600036, India;
sudheer@iitm.ac.in (K.P.S.); indunambi@iitm.ac.in (I.M.N.)
2 Agrosphere Institute (IBG-3), Forschungszentrum Jülich GmbH, 52425 Jülich, Germany;j.vanderborght@fz-juelich.de (J.V.); m.herbst@fz-juelich.de (M.H.)
3 Department of Environmental Sciences, University of California, Riverside, CA 92521, USA;jiri.simunek@ucr.edu
4 Department of Agricultural and Biological Engineering, Purdue University, West Lafayette, IN 47906, USA* Correspondence: beegumsahila@gmail.com; Tel.: +91-9497-308-126
1

Received: 26 February 2020; Accepted: 31 March 2020; Published: 2 April 2020

Abstract: Simulation models that describe the flow and transport processes of pesticides in soil andgroundwater are important tools to analyze how surface pesticide applications influence
groundwater quality. The aim of this study is to investigate whether the slow decline and the stablespatial pattern of atrazine concentrations after its ban, which were observed in a long-termmonitoring study of pesticide concentrations in the Zwischenscholle aquifer (Germany), could be
explained by such model simulations. Model simulations were carried out using MODFLOW model
coupled with the HYDRUS-1D package and MT3DMS. The results indicate that the spatial
variability in the atrazine application rate and the volume of water entering and leaving the aquifer
through lateral boundaries produced variations in the spatial distribution of atrazine in the aquifer.
The simulated and observed water table levels and the average annual atrazine concentrations were
found to be comparable. The long-term analysis of the simulated impact of atrazine applications inthe study area shows that atrazine persisted in groundwater even 20 years after its ban at an average
atrazine concentration of 0.035 µg/L. These results corroborate the findings of the previous
monitoring studies.
Keywords: MODFLOW with the HYDRUS-1D package and MT3DMS; pesticide contamination;
atrazine; groundwater contamination

1. Introduction
Pesticides are widely used all over the world to protect crops in agricultural fields. Uncontrolled
and improper pesticide applications have led to their presence in the environment exceeding
permissible limits [1–4]. It has been observed that groundwater is most affected by pesticide
applications in agricultural fields [5–9]. The pesticide contamination of groundwater is generally
assessed through pesticide monitoring programs (e.g., [10–12]). Once the information about pesticideconcentrations in the groundwater is obtained, it is crucial to analyze the spatiotemporal relations
between surface applications of pesticides and their groundwater concentrations. Such an analysis isimportant for several reasons: (a) to assess the role of non-point versus point scale sources ofcontamination, (b) to interpret observations of monitoring studies, (c) to identify vulnerable locationsin the aquifer, and (d) to evaluate the effectiveness of various regulatory measures for protecting
groundwater quality. The approaches that have been used for such assessments can be classified into
Water 2020, 12, 1019; doi:10.3390/w12041019

www.mdpi.com/journal/water

Water 2020, 12, 1019

2 of 24(a) geostatistical approaches [13], and the use of (b) spatially distributed one-dimensional soil models[14], (c) groundwater models [15], and (d) coupled soil and groundwater models [16,17]. The latter
type of modeling approaches can be used to assess the impact of various processes including
groundwater flow and solute transport, dilution, and decay, as well as the impact of the spatial
distribution of pesticide loadings to groundwater on the spatiotemporal patterns of pesticideconcentrations in the groundwater [18,19]. In this study, we investigate whether the slow decline andstable spatial patterns of atrazine concentrations after its ban, that were observed in a long-termmonitoring study of pesticide concentrations in the Zwischenscholle aquifer in Germany [20,21],
could be explained by such model simulations.
Atrazine was one of the most widely applied pesticides for weed control worldwide [20–24].
There are many cases of atrazine contamination of the soil and groundwater reported in different
parts of the world [25–28]. In Germany, atrazine was applied mainly in agricultural fields for maize
cultivation. This pesticide was banned in March 1991, since concentrations in groundwater anddrinking water were found to be exceeding the permissible limits. In Europe, according to theEuropean Council Directive 98/83/EC “The quality of water intended for human consumption” (1998)[29], the threshold limit of pesticide concentrations is 0.1 µg/L for a single compound and 0.5 µg/L
for the sum of all pesticides. Despite its ban, atrazine and its metabolite deethylatrazine wereobserved in German aquifers several years after its ban [24]. The Zwischenscholle aquifer in Germany
is one of many aquifers in the world contaminated with atrazine [22,23]. Due to a shallow
groundwater level and intensive agricultural land use in the Zwischenscholle aquifer, thevulnerability of this aquifer to pesticide contamination was found to be very high [24,30,31].
Vonberg et al. [20,21] conducted a detailed investigation of atrazine concentrations in the soil
and groundwater of the Zwischenscholle aquifer. In their study, the analysis of the observations froma monitoring campaign of atrazine concentrations from 1991 to 2011 showed that the groundwater
concentration of atrazine is still close to the threshold limit (0.1 µg/L even 20 years after its ban in
1991. The study demonstrated a considerable variation of atrazine concentrations in space and also
concluded that the spatial pattern of the concentration distribution was relatively stable during themonitoring period.
Centered on the concluding remarks made by Vonberg et al. [20,21], the main objective of thecurrent study was to investigate the atrazine contamination in the Zwischenscholle aquifer by
analyzing the persistent nature of the atrazine concentration and its spatiotemporal pattern using asimulation model. There exist various approaches and methods for analyzing the aquifer
contamination ranging from simple analytical to complex numerical methods. Analytical methods
facilitate a basic and rapid preliminary analysis of the groundwater contamination but are generally
associated with a number of simplifying assumptions regarding the groundwater system. They are
limited to simple aquifer geometries and homogeneous aquifer properties [32,33]. Some of themethods of groundwater contamination assessment are based on semi-quantitative approaches that
estimate the groundwater vulnerability by assigning weights to relevant parameters (hydrogeology,
topography, land use, soil, contaminant source, etc.) that control the movement of contaminants
through unsaturated and saturated soil zone (DRASTIC [34], SI [35], GOD [36], SGVI [37]). Some ofthese methods utilize remote sensing and GIS techniques (SINTACS [38]). There also exist coupledprocess-based and numerical approaches for groundwater contamination studies [39]. The stochastic
geostatistical data assimilation approach utilizes the solute concentration data from observation wells
to estimate the groundwater contamination and the plume distribution [40,41].
Simulation modeling of the groundwater system is another approach that can be used to analyze
the movement of water and contaminants in soil and groundwater by solving the fundamental
governing equations. Numerical methods can analyze irregular geometries and non-homogenous
structures of the aquifer. Simulation models that can analyze the movement of water and solute
through both unsaturated and saturated soil zones are important for a contaminant transport
analysis. There exist several numerical models that simulate these two zones independently [42–47]
or simultaneously [48–54]. The MODFLOW model with the HYDRUS-1D package and MT3DMS[53,54] was chosen for this study mainly due to its capability of combined modeling of water flow

Water 2020, 12, 1019

3 of 24and solute transport in the unsaturated and saturated soil zones and its well-adapted computational
complexity. The intention of using this simulation model in this study was to incorporate the effectof the highly heterogeneous characteristics of the study domain, varying land use, and several othercontrolling factors (e.g., locations and rates of atrazine applications, precipitation, potential
evapotranspiration, and changing water table elevations) on the transport of atrazine in the study
area, followed by simulating the long-term behavior of atrazine in the aquifer.
2. Materials and Methods
2.1. Study Area: Zwischenscholle Aquifer
The study area, the ‘Zwischenscholle aquifer,’ is located near Jülich in the Lower Rhine
Embayment, Germany. The Zwischenscholle aquifer is surrounded by Rurscholle separated by the
Rursprung fault in the southwest and by Erftscholle separated by the Rurrand fault in the northeast.
The Zwischenscholle area is lower than Erftscholle and higher than Rurscholle. Figure 1 shows thestudy area with geological faults and observation wells. Forest and agriculture are the major land
uses in this area. The western side is characterized by deciduous forest, whereas the eastern part isdominated by agricultural land use [55]. The unconfined or semi-confined aquifer consists of
Quaternary Rhine Maas sand and gravel sediments with a Reuver clay aquifer base. The groundwaterflow direction is from southeast to northwest (a blue arrow in Figure 1) with a mean hydraulicgradient between 0.1% and 0.2% [56]. The thickness of the aquifer varies from a few meters in thesouthwest to 35 m in the northeast. Land use, soil type, and soil hydraulic properties in the study
area are discussed in detail in Sections 2.3.2 and 2.3.3. Hydrogeological cross-sections A-A’ and B-B’of the aquifer are shown in Figure 2. It can be observed that the base of the upper aquifer is closer tothe surface in the southwest (a cross-section A-A’). There is an increase in the aquifer thickness in thenorth-western direction (a cross-section B-B’). The thin green line (the bottom of the Waal sedimentsin this area) represents the aquifer base considered in the current study. This is obtained from theinterpolation of the borehole data [55]. At some locations (e.g., southwest of cross-section A-A'),
where the Waal sediment layers were not present, the base of the aquifer is limited to depth with alower range of hydraulic conductivity. The green surface represents the Reuver clay, which is theaquifer base of a second (deeper) groundwater body. Four sides of the study area are represented in
Figure 1 as Sides 1, 2, 3, and 4.

Water 2020, 12, 1019

4 of 24

Figure 1. Study area with geological faults and observation wells.

Figure 2. Cross-sections A-A’ and B-B’ of the aquifer presented in Figure 1.

Water 2020, 12, 1019

5 of 24

2.2. MODFLOW with the HYDRUS-1D Package and MT3DMSIn MODFLOW with the HYDRUS-1D package [49,57] and MT3DMS [43], the one-dimensional
HYDRUS-1D [46] model for water flow and solute transport modeling in the unsaturated zone islinked to MODFLOW [58], a three-dimensional model for water flow modeling of the saturated zone,
and loosely coupled to MT3DMS (for solute transport modeling of the saturated zone) [53,54].
The HYDRUS-1D model [59] solves water flow in the unsaturated zone using the modified onedimensional Richards equation:

=

(ℎ)

ℎ

+1

− (ℎ)

(1)

where θ is the volumetric water content (dimensionless), h is the soil water pressure head (L), t is
time (T), z is the vertical coordinate (L) (positive upward), S is the sink term (T−1), and K(h) is theunsaturated hydraulic conductivity (LT−1).
HYDRUS-1D additionally simulates solute transport in variably saturated porous media using
the standard advection-dispersion transport equation of the form:

+

=

−

−∅

(2)

where c is the solution concentration (ML−3), s is the sorbed concentration (MM−1), D is the dispersion
coefficient (L2T−1), b is the bulk density of the porous medium (ML−3), q is the volumetric flux density(LT−1), which is obtained using the Darcy–Buckingham law, and ∅ is a sink–source term accounting
for various zero- and first-order or other reactions (ML−3T−1) [60]. In the current study, atrazine isconsidered to undergo the first-order degradation. Therefore,

∅= (

+

)

(3)

where λ is the first-order degradation rate (day−1). Root solute uptake is not considered in this study.
In MODFLOW, the three-dimensional movement of groundwater of a constant density is
described by the following partial differential equation:

ℎ

+

ℎ

ℎ

+

+

=

ℎ

(4)

where Kxx, Kyy, and Kzz are values of the hydraulic conductivity along the x, y, and z coordinate axes,
which are assumed to be parallel to the major axes of hydraulic conductivity (LT−1), h is thepotentiometric head (L), W is a volumetric flux per unit volume representing sources and/or sinks ofwater (T−1), SS is the specific storage of the porous material (L−1), and t is time (T).
MT3DMS simulate advection, dispersion/diffusion, and chemical reactions of contaminants inthe saturated zone using the following partial differential equation:

(

)

=

−

(

)+

+

(5)

where c is the dissolved solute concentration (ML−3), θis the porosity of the subsurface medium(dimensionless), t is time (T), xi is the distance along the respective Cartesian coordinate axis (L), Dij
is the hydrodynamic dispersion coefficient tensor (L2T−1), vi is the pore water velocity (LT−1), qs is thevolumetric flow rate per unit volume of the aquifer representing fluid sources (positive) and sinks(negative) (T−1), cs is the concentration of the source or sink flux (ML−3), and ∑is the chemical
reaction term (ML−3T−1).
In MODFLOW, the groundwater domain is discretized into finite difference grids. These grids
are then combined into multiple zones based on the similarities in soil characteristics, topographicalcharacteristics, boundary conditions, and the depth to groundwater. Each of these zones is thenassigned one soil profile that extends from the soil surface down to a depth, which should be below
the deepest possible water table level that can occur during the simulation (hereafter referred to asWater 2020, 12, 1019

6 of 24the ‘HYDRUS-1D profile’). The entire period of simulation in MODFLOW is divided into stressperiods, during which the input data for all external stresses (e.g., rainfall, potential
evapotranspiration, pumping, recharge, etc.) are constant. These stress periods are further divided
into time steps. Time steps are intervals in each stress period, during which the change in the aquifer
characteristics (e.g., groundwater head/water table level, groundwater flow, etc.) due to external
stresses is analyzed using the model. Water flow and solute transport are simulated in each
HYDRUS-1D profile for a duration equal to one MODFLOW time step by considering the averagewater table level from MODFLOW as the bottom pressure head boundary condition. The timeaveraged bottom flux from the HYDRUS-1D profile is given as the recharge flux to MODFLOW,
which then simulates groundwater flow. The average water table depth at the end of the MODFLOW
time step is assigned as the bottom pressure head boundary condition in the HYDRUS-1D profile forthe next MODFLOW time step [49,53,54,57]. The solute flux reaching the water table is further used
for solute transport modeling in the saturated soil zone using MT3DMS.
2.3. Data Available for Modeling
Aquifer data available for modeling and model validation include (i) information about
precipitation, potential evapotranspiration, boundary conditions, land use, and approximate
amounts of atrazine applied in agricultural fields from 1984 to 1993, (ii) atrazine concentrations indifferent observation wells in 1991, 1992, and from 2000 to 2011, except for the year 2004 [20,21], and(iii) water table fluctuations in selected observation wells in the study area from 1984 to 1993. The
information about land use and approximate quantities of atrazine applied in maize fields was
obtained using a questionnaire survey from the farmers [61]. Exact locations of maize fields andquantities of atrazine applied over the years in the study area are unknown since this information
was not recorded [61]. Therefore, the modeling is based on scenarios of atrazine applications that are
realistic in a statistical sense.
2.3.1. Climate Data
The climate data were obtained from the meteorological tower installed in the Jülich ResearchCentre [30]. Precipitation and potential evapotranspiration, calculated using the Penman–Monteithequation, from 1984 to 1993 are shown in Figure 3.

Figure 3. The daily precipitation and potential reference evapotranspiration data for the study areafrom 1984 to 1993.

Water 2020, 12, 1019

7 of 24

2.3.2. Land Use and Information about Atrazine Applications
The primary land use is agriculture. Since a complete map of land use was not available, astochastic distribution approach was used to generate possible hypothetical land use [30]. Figure 4
shows the spatial distribution of different land uses from 1984 to 1993. Winter wheat, sugar beet,pasture, maize, and potatoes represent about 50%, 35%, 8%, 4%, and 3%, respectively, of the total
agricultural land. The total area of maize fields is significantly smaller compared to other crops andis highly spatially spread in the study area.

Figure 4. Land use from 1984 to 1993.

The average application rate of atrazine in the agricultural field cultivated with maize is 0.0666
g/m [61]. Atrazine was applied once a year around May. Figure 5 shows the total mass of atrazine
applied in the study area and the area of the atrazine application from 1984 to 1993 [61].
2

(a)

(b)

Figure 5. The total mass of applied atrazine (a) and the area of the atrazine application in the study
area (b) from 1984 to 1993.

2.3.3. Soil Types and Soil Hydraulic Properties
The study area is characterized by loess sediments. Details regarding different soil types and soil
layers were available for the study area from the Forschungszentrum Jülich GmbH (as a part of the
MOSYRUR project [61], a project analyzing the water balance in the Rur basin). Soils in the study
area are classified into four soil types (silt loam, sandy loam, loam 1, and loam 2) and the parametersof the van Genuchten–Mualem soil hydraulic functions [62] were derived for different layers of theWater 2020, 12, 1019

8 of 24

four soil types. Table 1 shows the details about various soil types, soil layers, soil hydraulic properties(van Genuchten–Mualem analytical model parameters), and aquifer properties.
Table 1. Different types of soil, soil layers, soil hydraulic properties (the van Genuchten–Mualem
analytical model parameters) [62] and aquifer properties [30,61].

Soil
Type

Silt
loam
Sandy
loam
Loam 1
Loam 2Aquifer

LayersThickness(m)

Residual
Water
Content,
θr (-)

1
2
3
1
2
1
2
1
2
-

0.23
0.28
1.35
0.25
0.30
0.06
1.00
0.20
1.00
Variable

0.045
0.041
0.074
0.056
0.053
0.048
0.082
0.044
0.069
0.001

Saturated
Water
Content,θs (-)
0.463
0.392
0.359
0.382
0.314
0.461
0.358
0.395
0.298
0.184

Parameter
α in theSoil Water
Retention
Function,
(m−1)
2.16
1.60
1.42
11.2
9.33
2.10
4.32
2.04
4.16
10.0

Parametern in theSoil Water
Retention
Function,
(-)
1.35
1.35
1.30
1.35
1.40
1.34
1.33
1.34
1.33
2.10

Saturated
Hydraulic
Conductivity,
Ks (m day−1)

Tortuosity
Parameter l intheConductivity
Function (-)

0.080
0.032
0.017
1.38
0.891
0.074
0.067
0.049
0.110
134.

0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5

2.3.4. Simulation of Crop Transpiration and Root Water Uptake
The plant module SUCROS [63] was used by Herbst et al. [30] to estimate temporal changes in
the leaf area index (LAI) of various crops in the study area. The crop coefficient (Kc) is used to derive
specific crop evapotranspiration from potential reference evapotranspiration using the equation:
ETcrop = ETpot KC, where ETpot (LT−1) is potential evapotranspiration for the reference grass cover, ETcrop(LT−1) is plant-specific potential evapotranspiration. Temporal changes in the crop Kc factor followed
the Doorenbos and Pruitt [64] approach [61]. Actual transpiration from the soil is estimated in the
HYDRUS-1D package based on root water uptake characteristics. In HYDRUS-1D, root water uptake
is represented as extraction or a sink term distributed over the root zone. Root water uptake isevaluated using potential transpiration and the Feddes stress response function [65]. Parameters of
the Feddes stress response function used for different crops are given in Table 2. Rooting depths androot densities of various crops grown in the study area are given in Table 3.
Table 2. The Feddes stress response function parameters and maximum rooting depths of different
crops (P0 is the pressure head, below which roots start extracting water from the soil, POpt is thepressure head, below which roots extract water at the maximum possible rate, P2H is the value of thelimiting pressure head, below which roots can no longer extract water at the maximum rate of 0.5
cm/day, P2L is the value of the limiting pressure head, below which roots can no longer extract waterat the maximum rate of 0.1 cm/day, and P3 is the pressure head, below which root water uptake
ceases).
Crop
Winter wheat
Maize
PotatoSugar beet
Grassland
Forest

P0
(m)
0.000
−0.150−0.100−0.100−0.0011000.001

POpt(m)−0.010−0.300−0.250−0.250−0.025
1000

P2H
(m)
−5.00
−3.25
−3.20
−3.20
−2.00
−4.25

P2L
(m)
−9
−6
−6
−6
−8
−8

P3 (m)

Maximum Rooting Depth (m)

−160
−80
−160
−160
−80
−160

0.95
1.15
0.74
1.15
0.75
1.50

Water 2020, 12, 1019

9 of 24

Table 3. Relative rooting depths and root densities of different crops [55].
Relative
Rooting
Depth (-)
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0

Root Density; Value Characterizing the Depth Distribution of Root Water Uptake (m−1)
Winter Wheat

Maize

PotatoSugar Beet

Grassland

Forest

0.26
0.16
0.10
0.06
0.04
0.02
0.01
0.01
0.00
0.00
0.00

1.05
0.76
0.55
0.39
0.28
0.20
0.15
0.10
0.07
0.05
0.04

1.00
0.94
0.89
0.84
0.78
0.73
0.68
0.62
0.57
0.52
0.47

0.97
0.92
0.86
0.80
0.74
0.69
0.63
0.57
0.51
0.46
0.40

1.05
0.76
0.55
0.39
0.28
0.20
0.15
0.10
0.07
0.05
0.04

0.30
0.60
0.90
0.95
0.90
0.80
0.70
0.60
0.50
0.40
0.30

2.3.5. Solute Transport Parameters
Solute transport is described using the advection-dispersion equation. A linear sorption
isotherm is used to describe the sorption of atrazine. The sorption constant Kd (m3/kg), the first-orderdegradation rate, λ (day−1), the bulk density (kg/m3), and the longitudinal αL (m) and transverse αT(m) dispersivities for different soil types and soil depths are given in Table 4. It can be noted that thefirst-order degradation rate in the aquifer is very low (3.47 × 10-10 day−1) [55], which could be the main
reason for the long-term stable atrazine concentration in the groundwater [20,21].
Table 4. Solute transport parameters specific to soil and solute [55].
Soil
Type
Silt
loam
Sandy
loam
Loam 1
Loam 2Aquifer

LayersPorosity
(-)

1
2
3
1
2
1
2
1
2
-

0.509
0.434
0.434
0.509
0.434
0.510
0.440
0.509
0.434
0.184

Longitudinal
Dispersivity,
αL (m)
0.045
0.045
0.045
0.045
0.045
0.045
0.045
0.045
0.045
15.4

Transverse
Dispersivity,
αT (m)
0.00450.00450.00450.00450.00450.00450.00450.00450.0045
1.54

Bulk
Density,(kg/m3)
1350
1375
1375
1350
1375
1350
1375
1350
1375
1400

Adsorption
Coefficient,
Kd (m3/kg)
1.42 × 10-03
3.26 × 10-04
1.63 × 10-04
1.42 × 10-03
3.26 × 10-04
6.91 × 10-03
6.52 × 10-04
1.42 × 10-03
3.26 × 10-04
3.55 × 10-04

First-order
Degradation
Rate, λ (day−1)
3.47 × 10-03
1.73 × 10-03
1.04 × 10-03
3.47 × 10-03
1.73 × 10-03
3.47 × 10-03
1.73 × 10-03
3.47 × 10-03
1.04 × 10-03
3.47 × 10-10

2.4. Model Setup
2.4.1. Spatial and Temporal Discretizations
The study area of 21 km2 was divided into 504 MODFLOW grids each of a dimension of 200 ×
200 m (14 columns and 36 rows). Figure 6 shows the initial water table elevation, MODFLOW grids,
boundary conditions, the groundwater flow direction, and the top and bottom elevations of theaquifer at four boundary points. Considering the highly variable soil hydraulic properties, surface
elevations, water table elevations, and atmospheric boundary conditions in the study area, each
MODFLOW grid cell was assigned a single HYDRUS-1D profile. The top elevation of the HYDRUS1D profile was considered as the surface elevation corresponding to the MODFLOW grid cell
assigned to the profile. The bottom elevation of each HYDRUS-1D profile was set 2 m below the initialwater table elevation in the grid cell. This was further adjusted by performing trial runs that checked
for the possibility of the water table falling below the bottom of the HYDRUS-1D profile. EachHYDRUS-1D profile was then divided into finite elements based on the soil hydraulic properties andWater 2020, 12, 1019

10 of 24

the depth of the HYDRUS-1D profile. The ratio of the sizes of neighboring finite elements was always
less than 1.5 to ensure the convergence of the solution. The total number of finite elements varied
from 150 to 250 in all profiles. The number of stress periods and time steps used in this study is 14,609,
with the duration of each time step equal to 1 day. The functioning of the model is explained in
Section 2.2.

Figure 6. The initial water table elevation, MODFLOW grids, boundary condition, the generalgroundwater flow direction (the blue arrow), and the thickness of the aquifer at four boundary points.

2.4.2. Initial and Boundary Conditions
The water table elevation in the year 1984 (the beginning of the simulation) obtained byinterpolating groundwater table levels in the observation wells in the study area [61] was given as
the initial water table elevation (Figure 6). The atmospheric boundary condition was given at thesurface using precipitation and potential evapotranspiration fluxes. A no-flow boundary condition
was assumed at the bottom of the aquifer. Ten layers were considered for groundwater modelingusing MODFLOW and solute transport modeling using MT3DMS. The monthly varying water tablelevels were given as the variable head boundary conditions on all lateral boundaries of the study
domain. This data was obtained using a spatiotemporal interpolation procedure for defining theheads at the boundary at a 30-day time step from the measured groundwater levels in the aquifer[30].
2.4.3. Model Settings for Analyzing the Long-Term Atrazine Concentrations in Groundwater
To simulate the long-term atrazine concentrations in groundwater and also to validate themodel, the simulation period was extended until 2023. All information (precipitation, potential
evapotranspiration, boundary conditions, and land use), available for the first 10 years (1984–1993),
was repeated three times (from 1994 to 2023) during this exercise. This approach was used because
the data needed for the boundary conditions (the spatially and temporally interpolated groundwatertable levels) and land use were not available for the period after 1993. A similar approach was
adopted by Rakowski et al. [66] for developing future scenarios in groundwater modeling. The
information regarding precipitation and potential evapotranspiration from 1994 to 2023 could havebeen obtained either from the weather station in the study area or using weather forecasting.
However, since water table levels that could be used as boundary conditions were not available forthis period, and also since these water table levels depend on precipitation, it was decided to repeat
the information. During the 40-year simulation period, atrazine was applied during the first 20 years.
The evolution of atrazine concentrations during the first 20 years represents the change of atrazine
concentrations over time after the use of a substance in a certain region started. The last 20 years of
the simulation represent how the concentrations decrease after a ban on atrazine.

Water 2020, 12, 1019

11 of 24

2.4.4. Looping Boundary Condition
Since there is no specific information available regarding the inflow of solute with groundwater
across the lateral boundaries of the study domain, a methodology was developed in this study torepresent the effects of the entire basin on the study area (a part of the basin). Instead of considering
the inflow of clean water (water without atrazine) through the lateral aquifer boundaries, soluteconcentrations from MODFLOW cells in the layers which are saturated at the outflow side of thedomain were used in the corresponding active cells at the inflow side of the domain. This is further
referred to as a ‘looping boundary condition.’ In this case study, the looping boundary condition wasconsidered only in the main direction of groundwater flow, which is from south-east to north-west(from Side 4 to Side 2 in Figure 1). This boundary condition assumes that the average atrazine
concentration in water that percolates from the saturated soil zone to the downstream aquifer in the
study domain is equal to the average concentration of atrazine that reaches the groundwater in the
upstream part of the aquifer. Note that the use of this looping concentration boundary condition does
not prevent the atrazine mass in the aquifer to decrease over time. Due to incoming recharge,
outgoing water and atrazine fluxes are larger than incoming fluxes, and thus the total atrazine massin the aquifer decreases over time when recharge does not any longer contain atrazine (due to itsban). The source code of MT3DMS (http://hydro.geo.ua.edu/mt3d) was modified to incorporate thisspecial boundary condition.
3. Results and Discussion
3.1. Water Flow Modeling Using MODFLOW with the HYDRUS-1D Package
Figure 1 shows several observation wells in the study area represented by well identification
numbers. These are the wells, in which water elevations and atrazine concentrations were monitored[30]. Since there are only a few wells with information regarding the water table elevation with a hightemporal resolution, only four observation wells highlighted in red color (well numbers 20232, 20111,
20244, and 20255 in Figure 1) were used for model validation. The water table elevation data isavailable from October 1984 to November 1993, December 1983 to October 1986, December 1983 toNovember 1993, and December 1983 to November 1993 for wells 20232, 20111, 20244, and 20255,
respectively.
Figure 7 shows the comparison of model-simulated and observed water table elevations in the
study area. In general, observed and simulated water table elevations are found to be similar in all
four observation wells. Minor differences between observed and simulated water table elevations
could be the result of comparing observed point data with the average simulated water table positionin the grid cell (the grid cell area of 40,000 m2 in this study). Another plausible reason for the anomaly
between simulated and observed water tables could be the fact that the boundary condition provided
to the model was generated by spatiotemporal interpolation [30], which may also exhibit a certain
level of inaccuracy. Furthermore, the aquifer properties (hydraulic conductivity, specific yield, andporosity) may vary spatially. Additionally, a constant value of the specific yield was used in the
groundwater model simulations in the developed model. The specific yield can be considered
constant only when the aquifer response is linear, i.e., when the volume of water added or released
is linearly proportional to the water table fluctuation. The specific yield varies based on the transient
nature of the water release from the unsaturated soil zone, soil hydraulic properties, the depth to thegroundwater table, or hysteresis [67,68], which are factors not taken into account in this study.

Water 2020, 12, 1019

12 of 24

Figure 7. Comparison of water table fluctuations in observation wells (observed) and those simulatedusing MODFLOW with the HYDRUS-1D package (simulated).

3.2. Atrazine Concentrations in Groundwater
Figure 8 shows the average annual atrazine concentration in the groundwater obtained using
MODFLOW with the HYDRUS-1D package and MT3DMS. The horizontal dotted line shows thepermissible limit of atrazine concentrations, and the vertical line indicates the period when atrazine
applications in the field stopped. The atrazine concentration appeared in the groundwater about one
year after the start of its use. The atrazine concentration in the aquifer increased during the period of
its application and started slowly decreasing once its application was stopped. The rate at which theaverage solute concentration is increasing/decreasing in the aquifer depends on the rate of recharge
from the unsaturated soil zone to the saturated soil zone, the rate of atrazine application, the locationof the maize fields, lateral groundwater flow into and out of the saturated soil zone, and the looping
boundary condition.

Figure 8. Average annual atrazine concentrations in the aquifer simulated using MODFLOW with the
HYDRUS-1D package and MT3DMS.

Water 2020, 12, 1019

13 of 24

Figure 9 shows total net recharge (in m3/year) and the net recharge rate (in mm/year) between
the unsaturated and saturated soil zones in each year. Recharge depends on atmospheric boundary
conditions (precipitation, potential evaporation), land use, and soil hydraulic properties. It can beobserved that the net recharge varies from year to year, with a maximum rate of about 61 mm/year
and a minimum rate of about 11 mm/year.

Figure 9. Total net recharge and the net recharge rate between the unsaturated and saturated soil
zones in each year.

Figure 10a shows the total inflow and outflow of water through the four lateral sides of theaquifer. The amounts of water leaving and entering through the four lateral sides (Side 1 to Side 4 inFigure 1) are different because of groundwater gradients and differences in flow areas of each of thesesides. The maximum flow area is on Side 3, followed by Sides 1, 4, and 2. Most water is flowing into
the domain through Side 4 (Side 4_in), and most water is leaving the aquifer through Side 3 (Side3_out). In the looping direction (Sides 2–Side 4), more water is entering the domain through Side 4
than leaving through Side 2. Figure 10b shows volumes of water entering and leaving the aquifer
through four lateral sides divided by corresponding flow areas (area through which flux goes in orout), which represents average water fluxes. It can be observed that the relative volume (for the flowarea) of water entering the aquifer is highest on Side 4 (Side 4_in) and lowest on Side 2, while therelative volume of water leaving the aquifer is highest on Side 2 (Side 2_out) and lowest on Side 4(Side 4_out). This indicates that the major flow direction is from (Side 4) south-east to (Side 2) northwest, and the same direction is chosen for the solute transport looping boundary, as explained in
Section 2.4.

Water 2020, 12, 1019

14 of 24(a)(b)
Figure 10. Volumes of water entering and leaving through sides of the domain (a), and water fluxes
entering and leaving the aquifer through sides of the domain (b), (Side 1_in to Side 4_in represent
flow entering through particular lateral sides, and Side 1_out to Side 4_out represent flow leaving
through particular lateral sides).

Figure 11a shows the total amount of atrazine entering the groundwater through various sources
and leaving though various sinks. Various sources that add solute into the system include recharge
from the unsaturated zone to the saturated zone and lateral inflow through Side 4 (due to the looping
boundary condition). Lateral outflow from the aquifer acts as a sink. There was an increase in the
total amount of atrazine entering the saturated soil zone from year 3 to year 6 and from year 10 toyear 15 (Figure 11). This is because of an increase in the area of the atrazine application and also due
to leaching of adsorbed atrazine. There was a decrease in the amount of atrazine reaching theWater 2020, 12, 1019

15 of 24groundwater from year 7 to year 9 and after year 16. This is because of a decrease in the quantity of
applied atrazine and the reduction in recharge (Figure 9) during these periods. There was again an
increase in the amount of atrazine entering groundwater from year 22 to year 30 and a decrease
afterwards. This is because of the looping boundary condition to represent the influx of atrazine fromupstream areas.

(a)

(b)
Figure 11. The total mass of atrazine that is entering the aquifer through various sources and leaving
the aquifer through various sinks (a), and the total mass of atrazine that is entering the aquifer throughlateral inflow and recharge (b).

Water 2020, 12, 1019

16 of 24

Figure 11b shows the total mass of atrazine that is entering the aquifer through lateral flux (Side4) and recharge flux. Large amounts of atrazine are entering the domain through recharge during thefirst 20 years (when atrazine was applied in the area). When comparing Figure 11a,b, it can beobserved that the amount of atrazine entering the aquifer through lateral sides is smaller than theamount of atrazine leaving the aquifer through lateral sides during the first 12 years. This is becauseof the looping boundary condition and can be explained by Figure 12. Figure 12 shows the amountof atrazine leaving and entering through lateral sides. From Figure 12, it can be observed that only
small amounts of atrazine are leaving through Side 2 (Side 2_out) during the first 12 years. Since thelooping boundary is considered from Side 2 to Side 4, the amount of atrazine entering the lateral side
(Side 4) is also small during this period. During later years (after 12 years) the amount of atrazine
leaving through Side 2 increases, and so is the inflow through Side 4 in Figure 12.

Figure 12. The amount of atrazine entering and leaving the aquifer through four lateral sides (Side4_in represents atrazine volumes entering through Side 4, Side 1_out to Side 4_out represent atrazine
leaving through corresponding lateral sides, and Total Sides_out represents the amount of soluteleaving though all lateral sides).

From Figure 12, it can be observed that the amount of atrazine leaving through Sides 1 and 3 ishigher than through Side 2 in most years. This is because higher volumes of water are leaving throughSides 1 and 3 and also due to the location of atrazine applications. It can also be observed that theamount of solute leaving through Side 3 is always higher than through Side 2 and that the amount ofsolute entering through Side 4 due to the looping boundary condition corresponds to the amount ofsolute leaving the aquifer through Side 2. However, the amounts entering via Side 4 are larger than
the amount that leaves Side 2 because of the larger water flow into Side 4 than out of Side 2 (seeFigure 10a).
Figure 13 shows the spatial pattern of the depth-averaged atrazine concentration in the aquifer
from 1 to 40 years. In the early years of the simulation, based on the quantity and the location ofatrazine applied in the study area, a varying spatial pattern and a localized appearance of the atrazine
concentration was simulated in the aquifer. During the first 20 years of the simulation, distinct
atrazine plumes that emerged below the fields where atrazine was applied could be discerned. These
plumes were transported slowly with groundwater flow. Since no atrazine degradation wasconsidered in the aquifer, the concentrations in these plumes decreased slowly due to dilution causedby solute dispersion and mixing with recharge water that did not contain atrazine. Over time,Water 2020, 12, 1019

17 of 24atrazine spread more and more in the aquifer, and it was still observed in the aquifer several years
after year 20 (the year after which atrazine applications stopped). The spatial pattern of the atrazine
distribution after year 20 remains quite stable over the next 20 years and is mainly influenced by thelooping boundary condition (Figure 13).

Figure 13. The spatial pattern of the atrazine concentration (µg/L) in the aquifer from year 1 to year
40 (a logarithmic color scale).

At the beginning of the simulation (year 1 to year 6), the solute concentration was increasing in
the aquifer (Figures 8 and 13) as more and more atrazine was applied at various locations in the study
area. During the next few years (year 7 to year 9), there was a slight decrease in the atrazine
concentration mainly because of a reduction in atrazine applications (Figure 5) and a reduction inrecharge during this period (Figure 9). This was also because of larger lateral outflow during thisperiod. It should also be noted that the looping boundary condition was considered only in the main
direction of groundwater flow, whereas lateral outflow was also observed in other directions during
this period, which carried solute out of the system and which resulted in a decrease in theconcentration in this period.
There was again an increase in the atrazine concentration from year 10 to year 16 (Figures 8 and13), mainly because of the combined effect of larger atrazine applications (Figure 5) and larger
recharge (Figure 9). Larger recharge has a higher potential of leaching into the groundwater moreatrazine either applied at the soil surface or already present in the unsaturated soil zone. Severeleaching occurs when a large amount of atrazine is applied at the soil surface simultaneously with

Water 2020, 12, 1019

18 of 24meteorological conditions that result in larger recharge to the saturated soil zone. A similar case exists
during years 13 to 15, resulting in higher groundwater concentrations (Figures 8 and 11a).
There is a reduction in atrazine concentrations from year 16 to year 17 because of lower recharge
and larger lateral outflow. The average atrazine concentration in the aquifer after its application was
stopped in the year 20 was decreasing at a very slow rate. This is similar to the observation of Vonberg
et al. [20,21]. The amount of atrazine in the aquifer after 20 years of its application is influenced byavailable atrazine in the soil and groundwater, subsequent leaching, and atmospheric boundary
conditions. However, the amount of atrazine that enters in the aquifer via recharge declines rapidly
after the application of atrazine in the area was stopped, suggesting that the stock of atrazine in thesoil that could contribute to leaching is depleted rapidly because of the degradation in the soil layers.
A small increase in the average atrazine concentration was observed from year 27 to year 32 becauseof the inflow of atrazine due to the looping boundary condition adopted in this case study.
Six observation wells were chosen to analyze the spatial and temporal evolution of the atrazine
concentration in the study area after the atrazine application was stopped. The locations of theobservation wells, along with the MODFLOW grids, are shown in Figure 14b. Figure 14a shows
atrazine concentrations in the observation wells after 20 years of its application. When examining
concentrations in the observation wells after the atrazine application was stopped (Figure 14a), it isobserved that the concentration in observation wells 2 and 3 is decreasing while it is increasing inobservation wells 1, 5, and 6. An increase in concentrations is caused by the movement of atrazine
plumes to these locations. The concentration in the observation well 4 remains more or less constant.
In general, when examining selected observation wells, it is observed that the temporal evolution ofatrazine concentrations is different at different locations in the study area. This is similar to theobservation made by Vonberg et al. [20,21], who concluded that the study area has a considerable
variation of atrazine concentrations in space.

(a)(b)
Figure 14. Atrazine concentrations in six observation wells after 20 years of atrazine application (a),and the locations of the observation wells (b).

3.3. Comparison of Observed and Model-Simulated Average Annual Atrazine Concentrations
Vonberg et al. [20] conducted a detailed study of atrazine concentrations in groundwater of the
Zwischenscholle aquifer based on observations in monitoring wells. This study showed that several
modifications and changes were adopted during the monitoring program. The number and locationof monitoring wells, the number of samples collected, and the number of detects and non-detects inthe samples collected were different over time. Observation wells with small or no atrazine
concentrations were excluded, and some wells with higher concentrations were added during themonitoring program [20]. Figure 15 shows the total number of wells, the number of samples, and thenumber of detects and non-detects. This data collected for the study by Vonberg et al. [20] was not
sufficient to make a conclusive statement about the change in the atrazine concentration in theWater 2020, 12, 1019

19 of 24

aquifer. In addition to this, there is only limited knowledge about the actual location, quantity, andtime of the atrazine application in the study area. Therefore, a direct comparison between simulated
and observed solute concentrations with respect to the location and time of atrazine applications is
not possible. Nonetheless, in order to validate the performance of the model, the average annual
atrazine concentration in groundwater simulated by the model is compared here with the averageannual atrazine concentration obtained from field measurements. The years from 2000 to 2011
correspond to 29th to 40th year in the simulation study. Observation data is not available for the year
2004 (33rd year in the simulation study).

Figure 15. The number of observation wells, samples collected, samples with detects and non-detectsfrom the monitoring program in the Zwischenscholle aquifer [20,21].

Figure 16 shows boxplots of the annual atrazine concentration observed in monitoring wells, andmodel-simulated concentrations along with the mean of observed and simulated concentrations
between years 29 and 40 (except for the year 33). The average annual atrazine concentration in theaquifer obtained using MODFLOW with the HYDRUS-1D package and MT3DMS over the entire
simulation period was found to be below the threshold level suggested by the European Council
Directive 98/83/EC “The quality of water intended for human consumption” (1998) [29] (0.1 µg/L). The
means of observed and simulated atrazine concentrations were found to be comparable. The meanatrazine concentration in observation wells is closer to the threshold level than the model-simulated
values, which is mainly because of the addition of monitoring wells, which showed higher
concentrations, and the exclusion of wells, which showed smaller concentrations in the monitoring
program as explained earlier in this section. Similar observations were made in several other studies on
the fate of atrazine in the Zwischenscholle aquifer and other parts of Germany, showing concentrations
close to the threshold limit [21,24,69] even after several years of atrazine ban in the country.

Water 2020, 12, 1019

20 of 24

Figure 16. Boxplots of the annual atrazine concentrations observed in monitoring wells (represented
in blue color) and simulated (represented in brown color), along with the arithmetic mean of observed
and simulated atrazine concentrations (a logarithmic scale).

The boxplots of the observed and simulated atrazine concentrations show that simulated ranges
are narrower than those observed, likely because of the mixing of local concentrations due todispersion and numerical averaging over 200 × 200 m grid cells. The groundwater sampling is
restricted to the volume of groundwater that is sampled in a single observation well. In groundwatertransport models, a macroscopic dispersion coefficient is used to represent the effect of aquifer
heterogeneity on the solute plume spreading. However, advection dominated spreading can be muchlarger than mixing and dilution controlled by local-scale dispersion. Therefore, local measured
concentrations may vary much more than simulated concentrations because local transport
heterogeneities are not averaged out by local measurements.
4. Concluding Remarks
The atrazine contamination in the Zwischenscholle aquifer was analyzed in this study using
MODFLOW with the HYDRUS-1D package and MT3DMS. The simulated and observed water tablelevels and atrazine concentrations were comparable. Model simulations indicate that the variability
in the atrazine application rate and the volume of water entering and leaving the aquifer throughvarious sinks and sources led to the variations in the amount of atrazine and its spread in the aquifer.
The impact of these different fluxes and their variability was assessed using the model. The longterm analysis of the impact of the atrazine application in the study area shows that the atrazine ispersisting in groundwater even after 20 years of its ban. Similar observations were made by Vonberg
et al. [20,21] after the detailed investigation of observed atrazine concentrations in the observation
wells in the study area.
Though the simulation model can, in general, model the effects of aquifer heterogeneity andvarious other factors that affect the movement of atrazine in the aquifer, the model was not able toaccurately explain the heterogeneous distribution of atrazine concentrations in the aquifer (observedin the monitoring study of Vonberg et al. [20,21]) because of the limited knowledge about the exact
quantity and location of atrazine applications in the study area.
The model could not resolve the small-scale transport variability in the aquifer sinceconcentrations and flow velocities were averaged over 200 × 200 m grid cells. Concentrations
monitored in groundwater monitoring wells are influenced by these small-scale heterogeneities. The
variation of concentrations observed in different monitoring wells can, therefore, be larger than thespatial variation of simulated concentrations. Simulations explained that an observed increase inatrazine concentrations in individual observation wells after the stop of atrazine applications were

Water 2020, 12, 1019

21 of 24

the result of transport processes in the aquifer, in combination with spatial variations of atrazineleaching from the soil to the aquifer that generated spatial heterogeneity of atrazine concentrationsin the groundwater. Presented simulations demonstrate the great modeling capabilities of thecoupled modeling system involving MODFLOW with the HYDRUS-1D package and MT3DMS.
Author Contributions: S.B., J.V., J.Š. and M.H. contributed to the conceptualization, methodology, analysis, andmanuscript preparation; J.V. and M.H. contributed to the resources and data, K.P.S. and I.M.N. contributed toreviewing and editing. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.
Acknowledgments: The first author S.B. would like to acknowledge the support from the Indo-German Centre
for Sustainability (IGCS) Scholarship for research exchanges in India and Germany for the research stay at
Agrosphere Institute (IBG-3), Forschungszentrum Jülich GmbH, 52425, Jülich, Germany for pursuing this work.
Conflicts of Interest: The authors declare no conflict of interest.

References
1.

2.

3.

4.
5.
6.

7.

8.
9.
10.
11.
12.

13.
14.
15.
16.

Huang, F.; Li, Z.; Zhang, C.; Habumugisha, T.; Liu, F.; Luo, X. Pesticides in the typical agricultural
groundwater in Songnen plain, northeast China: Occurrence, spatial distribution and health risks. Environ.
Geochem. Health 2019, 41, 2681–2695.
Ali, N.; Khan, S.; Ur Rahman, I.; Muhammad, S. Human health risk assessment through consumption oforganophosphate pesticide-contaminated water of Peshawar basin, Pakistan. Expo. Health 2018, 10, 259–
272.
Gupta, S.D.; Mukherjee, A.; Bhattacharya, J.; Bhattacharya, A. An overview of agricultural pollutants andorganic contaminants in groundwater of India. In Groundwater of South Asia; Springer: Singapore, 2018; pp.
247–255.
El Alfy, M.; Faraj, T. Spatial distribution and health risk assessment for groundwater contamination from
intensive pesticide use in arid areas. Environ. Geochem. Health 2017, 39, 231–253.
Chaza, C.; Sopheak, N.; Mariam, H.; David, D.; Baghdad, O.; Moomen, B. Assessment of pesticidecontamination in Akkar groundwater, northern Lebanon. Environ. Sci. Pollut. Res. 2018, 25, 14302–14312.
Pan, H.; Lei, H.; He, X.; Xi, B.; Han, Y.; Xu, Q. Levels and distributions of organochlorine pesticides in thesoil–groundwater system of vegetable planting area in Tianjin City, Northern China. Environ. Geochem.
Health 2017, 39, 417–429.
Sieczka, A.; Bujakowski, F.; Falkowski, T.; Koda, E. Morphogenesis of a floodplain as a criterion forassessing the susceptibility to water pollution in an agriculturally rich valley of a lowland river. Water 2018,10, 399.
Segerson, K. Liability for groundwater contamination from pesticides. Econ. Water Qual. 2018, 19, 227–243.
Srivastav, A.L. Chemical fertilizers and pesticides: Role in groundwater contamination. In Agrochemicals
Detection, Treatment and Remediation; Elsevier: Amsterdam, The Netherlands, 2020; pp. 143–159.
Worrall, F.; Besien, T. The vulnerability of groundwater to pesticide contamination estimated directly fromobservations of presence or absence in wells. J. Hydrol. 2005, 303, 92–107.
Thapinta, A.; Hudak, P.F. Use of geographic information systems for assessing groundwater pollution
potential by pesticides in Central Thailand. Environ. Int. 2003, 29, 87–93.
Papadopoulou-Mourkidou, E.; Karpouzas, D.G.; Patsias, J.; Kotopoulou, A.; Milothridou, A.; Kintzikoglou,
K.; Vlachou, P. The potential of pesticides to contaminate the groundwater resources of the Axios river
basin in Macedonia, Northern Greece. Part I. Monitoring study in the north part of the basin. Sci. Total
Environ. 2004, 321, 127–146.
Leterme, B.; Vanclooster, M.; Rounsevell, M.D.A.; Bogaert, P. Discriminating between point and non-point
sources of atrazine contamination of a sandy aquifer. Sci. Total Environ. 2006, 362, 124–142.
Tiktak, A.; De Nie, D.S.; Garcet, J.D.P.; Jones, A.O.; Vanclooster, M. Assessment of the pesticide leachingrisk at the Pan-European level. The EuroPEARL approach. J. Hydrol. 2004, 289, 222–238.
Loague, K.; Soutter, L.A. Desperately seeking a cause for hotspots in regional-scale groundwater plumes
resulting from non-point source pesticide applications. Vadose Zone J. 2006, 5, 204–221.
Franke, H.-J.; Teutsch, G. Stochastic simulation of the regional pesticide transport including theunsaturated and the saturated zone. Ecol. Model. 1994, 75, 529–539.

Water 2020, 12, 1019

17.

18.
19.
20.
21.

22.

23.
24.
25.
26.
27.

28.
29.
30.

31.

32.
33.
34.

35.
36.

37.
38.

22 of 24

Kupfersberger, H.; Klammler, G.; Schumann, A.; Brückner, L.; Kah, M. Modeling subsurface fate of smetolachlor and metolachlor ethane sulfonic acid in the Westliches Leibnitzer Feld aquifer. Vadose Zone J.
2018, 17, 170030, doi:10.2136/vzj2017.01.0030.
Honeycutt, R. Mechanisms of Pesticide Movement into Ground Water; CRC Press: Boca Raton, FL, USA, 2018;
ISBN 1351082795.
Jones, R.L. Modeling the degradation and movement of agricultural chemicals in ground water. In
Mechanisms of Pesticide Movement into Ground Water; CRC Press: Boca Raton, FL, USA, 2018; 117–124.
Vonberg, D.; Vanderborght, J.; Cremer, N.; Pütz, T.; Herbst, M.; Vereecken, H. 20 years of long-term
atrazine monitoring in a shallow aquifer in western Germany. Water Res. 2014, 50, 294–306.
Vonberg, D.; Hofmann, D.; Vanderborght, J.; Lelickens, A.; Köppchen, S.; Pütz, T.; Burauel, P.; Vereecken,
H. Atrazine soil core residue analysis from an agricultural field 21 years after its ban. J. Environ. Qual. 2014,
43, 1450.
Farlin, J.; Drouet, L.; Gallé, T.; Pittois, D.; Bayerle, M.; Braun, C.; Maloszewski, P.; Vanderborght, J.; Elsner,
M.; Kies, A. Delineating spring recharge areas in a fractured sandstone aquifer (Luxembourg) based onpesticide mass balance. Hydrogeol. J. 2013, 21, 799–812.
Gutierrez, A.; Baran, N. Long-term transfer of diffuse pollution at catchment scale: Respective roles of soil,
and the unsaturated and saturated zones (Brévilles, France). J. Hydrol. 2009, 369, 381–391.
Tappe, W.; Groeneweg, J.; Jantsch, B. Diffuse atrazine pollution in German aquifers. Biodegradation 2002,
13, 3–10.
Dunlap, K.; Brown, A. Atrazine Levels in Water, Sediment, and Amphibian Tissue Samples from Selected Ponds
in Westernmost Kentucky; Murray State University: Murray, KY, USA, 2018.
Guzzella, L.; Pozzoni, F.; Giuliano, G. Herbicide contamination of surficial groundwater in Northern Italy.
Environ. Pollut. 2006, 142, 344–353.
Jayachandran, K.; Steinheimer, T.R.; Somasundaram, L.; Moorman, T.B.; Kanwar, R.S.; Coats, J.R.
Occurrence of atrazine and degradates as contaminants of subsurface drainage and shallow groundwater.
J. Environ. Qual. 1994, 23, 311–319.
Wehtje, G.; Leavitt, J.R.C.; Spalding, R.F.; Mielke, L.N.; Schepers, J.S. Atrazine contamination ofgroundwater in the Platte Valley of Nebraska from non-point sources. Sci. Total Environ. 1981, 21, 47–51.
Regulation, E.U. Council Directive 98/83/EC of 3 November 1998 on the quality of water intended forhuman consumption. Off. J. Eur. Commun. 1998, 41, 32–54.
Herbst, M.; Hardelauf, H.; Harms, R.; Vanderborght, J.; Vereecken, H. Pesticide fate at regional scale:
Development of an integrated model approach and application. Phys. Chem. Earth Parts A/B/C 2005, 30, 542–
549.
Mouvet, C.; Baran, N.; Thiery, D.; Gutierrez, A.; Darsy, C.; Dubus, I.G.; Falkiewicz, W.; Ritsema, C.; Dekker,
L.; Normand, B. Pesticides in European groundwaters: Detailed study of representative aquifers andsimulation of possible evolution scenarios (PEGASE). In Proceedings of the EU Workshop: The Functioning
and Management of the Water-Soil-System at River-Basin Scale: Diffuse Polltuion and Point Sources,
Orléans, France, 26–28 November 2003; pp. 248–263.
Suffet, I.H.; Gibs, J.; Coyle, J.A.; Chrobak, R.S.; Yohe, T.L. Applying analytical techniques to solve
groundwater contamination problems. J. Am. Water Works Assoc. 1985, 77, 65–72.
Walton, W.C. Analytical Groundwater Modeling Flow and Contaminant Migration; Lewis Publishers: Chelsea,
MI, USA, 1989.
Aller, L. DRASTIC: A Standardized System for Evaluating Ground Water Pollution Potential Using Hydrogeologic
Settings; Robert, S., Ed.; Kerr Environmental Research Laboratory, Office of Research and Development,
US Environmental Protection Agency: Washington, DC, USA, 1985.
Bartzas, G.; Tinivella, F.; Medini, L.; Zaharaki, D.; Komnitsas, K. Assessment of groundwater contamination
risk in an agricultural area in north Italy. Inf. Process. Agric. 2015, 2, 109–129.
Djoudi, S.; Boulabiez, F.; Pistre, S.; Houha, B. Assessing Groundwater Vulnerability to Contamination in a
Semi-Arid Environment Using DRASTIC and GOD Models, Case of F’kirina Plain, North of Algeria. IOSR
J. Environ. Sci., Toxicol. Food Technol. 2019, 13, 39–44.
Al-Adamat, R.; Al-Shabeeb, A.A.-R. A simplified method for the assessment of groundwater vulnerability
to contamination. J. Water Resour. Protect. 2017, 9, 305.
Civita, M.; De Maio, M.; Ubertini, L. Valutazione e cartografia automatica della vulnerabilità degli acquiferi
all’inquinamento con il sistema parametrico Sintacs R5: A new parametric system for the assessment and automaticWater 2020, 12, 1019

39.

40.
41.

42.

43.

44.

45.

46.

47.
48.
49.
50.
51.

52.
53.
54.

55.

56.
57.

58.

59.

23 of 24

mapping of ground water vulnerability to contamination; Pitagora: Bologna, Italy, 2000, ISBN 8837112319.
Fusco, F.; Allocca, V.; Coda, S.; Cusano, D.; Tufano, R.; De Vita, P. Quantitative assessment of specificvulnerability to nitrate pollution of shallow alluvial aquifers by process-based and empirical approaches.
Water 2020, 12, 269.
Michalak, A.M.; Shlomi, S. A geostatistical data assimilation approach for estimating groundwater plume
distributions from multiple monitoring events. Geophys. Monogr. Am. Geophys. Union 2007, 171, 73.
McLean, M.I.; Evers, L.; Bowman, A.W.; Bonte, M.; Jones, W.R. Statistical modelling of groundwatercontamination monitoring data: A comparison of spatial and spatiotemporal methods. Sci. Total Environ.
2019, 652, 1339–1346.
Lappala, E.G.; Healy, R.W.; Weeks, E.P. Documentation of Computer Program VS2D to Solve the Equations of
Fluid Flow in Variably Saturated Porous Media; Department of the Interior, US Geological Survey: Reston, VA,
USA, 1987; p. 83.
Zheng, C.; Wang, P.P. MT3DMS: A Modular Three-Dimensional Multispecies Transport Model for Simulation of
Advection, Dispersion, and Chemical Reactions of Contaminants in Groundwater Systems; Documentation andUser’s Guide; Alabama University: Tuscaloosa, AL, USA, 1999.
Xu, T.; Sonnenthal, E.; Spycher, N.; Pruess, K. TOUGHREACT: A New Code of the TOUGH Family for NonIsothermal Multiphase Reactive Geochemical Transport in Variably Saturated Geologic Media; University ofCalifornia: Berkeley, CA, USA, 2003.
Healy, R.W. Simulation of Solute Transport Invariably Saturated Porous Media with Supplemental Information on
Modifications to the US Geological Survey’s Computer Program VS2D.; Department of the Interior, US
Geological Survey: Reston, VA, USA, 1990; p. 90.
Šimůnek, J.; van Genuchten, M.T.; Šejna, M. The HYDRUS-1D software package for simulating the onedimensional movement of water, heat, and multiple solutes in variably-saturated media. Univ. Calif.
Riverside Res. Rep. 2005, 3, 1–240.
Šimůnek, J.; Bradford, S.A. Vadose zone modeling: Introduction and importance. Vadose Zone J. 2008, 7,
581–586.
Refsgaard, J.; Storm, B. Computer models of watershed hydrology, MIKE SHE. Water Resour. Publ. 1995,
809–846.
Twarakavi, N.K.C.; Simunek, J.; Seo, S. Evaluating interactions between groundwater and vadose zone
using the HYDRUS-based flow package for MODFLOW. Vadose Zone J. 2008, 7, 757–768.
Brunner, P.; Simmons, C.T. HydroGeoSphere: A fully integrated, physically based hydrological model.
Groundwater 2012, 50, 170–176.
Kalbacher, T.; Delfs, J.O.; Shao, H.; Wang, W.; Walther, M.; Samaniego, L.; Schneider, C.; Kumar, R.;
Musolff, A.; Centler, F.; et al. The IWAS-ToolBox: Software coupling for an integrated water resourcesmanagement. Environ. Earth Sci. 2012, 65, 1367–1380.
Morway, E.D.; Niswonger, R.G.; Langevin, C.D.; Bailey, R.T.; Healy, R.W. Modeling variably saturated
subsurface solute transport with MODFLOW-UZF and MT3DMS. GroundWater 2013, 51, 237–251.
Beegum, S.; Šimůnek, J.; Szymkiewicz, A.; Sudheer, K.P.; Nambi, I.M. Implementation of solute transport
in the vadose zone into the “HYDRUS Package for MODFLOW.” Groundwater 2019, 57, 392–408.
Beegum, S.; Šimůnek, J.; Szymkiewicz, A.; Sudheer, K.P.; Nambi, I.M. Updating the coupling algorithm
between HYDRUS and MODFLOW in the HYDRUS Package for MODFLOW. Vadose Zone J. 2018, 17,180034.
Mouvet, C.; Albrechtsen, H.J.; Baran, N.; Chen, T.; Clausen, L.; Darsy, C. PEGASE. Pesticides in European
Groundwaters: Detailed Study of Representative Aquifers and Simulation of possible Evolution Scenarios; Final
Report of the European Project; UW Umweltwirtschaft GmbH: Stuttgart, Germany, 2004.
Lahmeyer Site appraisal of geology, hydrogeology and soil mechanics (Standortgutachten Geologie,
Hydrogeologie und Bodenmechanik). Unpublished work, 1984.
Seo, H.S.; Šimůnek, J.; Poeter, E.P. Documentation of the HYDRUS Package for MODFLOW -2000, the US
Geological Survey Modular Ground Water Model. GWMI 2007-01, IGWMC-International Ground Water ModelingCenter; Colorado School of Mines: Golden, CO,USA, 2007; p. 96.
Harbaugh, A.W.; Banta, E.R.; Hill, M.C.; Mcdonald, M.G. MODFLOW-2000, the US Geological Survey
modular ground-water model: User guide to modularization concepts and the ground-water flow process.
USA Geol. Surv. 2000, 0–92.
Šimůnek, J.; van Genuchten, M.T.; Šejna, M. Recent developments and applications of the HYDRUS

Water 2020, 12, 1019

60.
61.
62.
63.

64.
65.
66.
67.
68.
69.

24 of 24computer software packages. Vadose Zone J. 2016, 15, doi:10.2136/vzj2016.04.0033.
Šimůnek, J.; van Genuchten, M.T. Modeling nonequilibrium flow and transport processes using HYDRUS.
Vadose Zone J. 2008, 7, 782–797.
Bogena, H. MOSYRUR: Water Balance Analysis in the Rur Basin; Heye Bogena, Forschungszentrum,
Zentralbibliothek: Julich, Germany, 2005.
van Genuchten, M.T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils.
Soil Sci. Soc. Am. J. 1980, 44, 892–898.
Spitters, C.J.T.; Van Keulen, H.; Van Kraalingen, D.W.G. A simple and universal crop growth simulator:
SUCROS87. In Simulation and Systems Management in Crop Protection; Rabbinge, R., Ward, S.A., Laar, H.H.
van, Eds.; Pudoc: Wageninge, Netherlands, 1989; pp. 147–181. ISBN 9022008991.
Doorenbos, J.; Pruitt, W.O. Crop Water Requirements. FAO Irrigation and Drainage Paper 24; Land and Water
Development Division; FAO: Rome, Italy, 1977; p. 144.
Feddes, R.A.; Kowalik, P.J.; Zaradny, H. Simulation of Field Water Use and Crop Yield; John Wiley and Sons:
New York, NY, USA, 1978.
Rakowski, P.; Knowling, M.J. Heretaunga Aquifer Groundwater Model: Development Report; HBRC Report No
RM18-14-4997; Hawke's Bay Regional Council: Napier, New Zealand, 2018.
Nachabe, M.H. Analytical expressions for transient specific yield and shallow water table drainage. Water
Resour. Res. 2002, 38, 11.
Said, A.; Nachabe, M.; Ross, M.; Vomacka, J. Methodology for estimating specific yield in shallow water
environment using continuous soil moisture data. J. Irrig. Drain. Eng. 2005, 131, 533–538.
Jablonowski, N.D.; Köppchen, S.; Hofmann, D.; Schäffer, A.; Burauel, P. Persistence of 14C-labeled atrazine
and its residues in a field lysimeter soil after 22 years. Environ. Pollut. 2009, 157, 2126–2131.
© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open
access article distributed under the terms and conditions of the Creative
Commons
Attribution(CC
BY)
license(http://creativecommons.org/licenses/by/4.0/).