J. Hydrol. Hydromech., 66, 2, 246–256
DOI: 10.2478/johh-2018-0005

Simulations of freshwater lens recharge and salt/freshwater interfaces using
the HYDRUS and SWI2 packages for MODFLOW
Adam Szymkiewicz1*, Anna Gumuła-Kawęcka1, Jirka Šimůnek2, Bertrand Leterme3, Sahila Beegum2,
Beata Jaworska-Szulc1, Małgorzata Pruszkowska-Caceres1, Wioletta Gorczewska-Langner1,
Rafael Angulo-Jaramillo4, Diederik Jacques31Gdańsk University of Technology, Faculty of Civil and Environmental Engineering, ul. Narutowicza 11, 80-233 Gdańsk, Poland.
University of California Riverside, 900 University Ave., Riverside, CA 92521, USA.
3
Engineered and Geosystems Analysis, Institute for Environment, Health and Safety, Belgian Nuclear Research Centre, Boeretang 200,
2400 Mol, Belgium.
4
Laboratoire d'Ecologie des Hydrosystèmes Naturels et Anthropisés (LEHNA) UMR 5023, 3, rue Maurice Audin, 69518 Vaulx-en-Velin,France.
*Corresponding author. E-mail: adams@pg.gda.pl
2

Abstract: The paper presents an evaluation of the combined use of the HYDRUS and SWI2 packages for MODFLOW
as a potential tool for modeling recharge in coastal aquifers subject to saltwater intrusion. The HYDRUS package for
MODFLOW solves numerically the one-dimensional form of the Richards equation describing water flow in variablysaturated media. The code computes groundwater recharge to or capillary rise from the groundwater table whileconsidering weather, vegetation, and soil hydraulic property data. The SWI2 package represents in a simplified way
variable-density flow associated with saltwater intrusion in coastal aquifers. Combining these two packages within the
MODFLOW framework provides a more accurate description of vadose zone processes in subsurface systems withshallow aquifers, which strongly depend upon infiltration. The two packages were applied to a two-dimensional problem
of recharge of a freshwater lens in a sandy peninsula, which is a typical geomorphologic form along the Baltic and theNorth Sea coasts, among other places. Results highlighted the sensitivity of calculated recharge rates to the temporal
resolution of weather data. Using daily values of precipitation and potential evapotranspiration produced average
recharge rates more than 20% larger than those obtained with weekly or monthly averaged weather data, leading todifferent trends in the evolution of freshwater-saltwater interfaces. Root water uptake significantly influenced both therecharge rate and the position of the freshwater-saltwater interface. The results were less sensitive to changes in soil
hydraulic parameters, which in our study were found to affect average yearly recharge rates by up to 13%.
Keywords: Groundwater recharge; Freshwater lenses; HYDRUS; SWI2; MODFLOW; Vadose zone.
INTRODUCTION
Groundwater recharge through the vadose zone is an important factor affecting salinity patterns in coastal aquifers. In
particular, in view of possible climate change, increased water
demands and sea level rise, the development of freshwater
lenses supplied with infiltrated rainwater on top of saline aquifers is a subject of intense research (e.g., Werner et al., 2013).
A large number of contributions have focused on freshwater
lenses in deltaic regions, as well as on islands, peninsulas, or
sandbars located in tropical, arid, and temperate climatic zones(e.g., Comte et al., 2014; Contractor and Jenson, 2000; Changet al., 2016; De Louw et al., 2011; Eeman et al., 2017; Holding
and Allen, 2015; Houben and Post, 2016; Illangasekare et al.,2006; Jocson et al., 2002; Mahmoodzadeh et al., 2014; Sadurski
et al., 1987; Sinclair et al., 2016; Stuyfzand, 2016; Sulzbacher
et al., 2012).
The actual amount of recharge reaching the groundwater table depends on several factors, including time-variable weather
conditions, land use, hydraulic properties of the vadose zone,
topographic relief, and the groundwater table depth. Thus,
recharge is inherently variable in space and time and this variability may have a significant influence on freshwater-saltwater
dynamics in shallow coastal aquifers (Mollema and Antonelli,
2013; Trglavcnik et al., 2016; Vandenbohede et al., 2014).
Other studies highlighted the importance of such factors as
water uptake by tree roots (Comte et al., 2014) or hydraulic246

properties and water storage in the vadose zone (Contractor andJenson, 2000; Jocson et al., 2002). There is thus a need to include seasonal fluctuations of recharge in groundwater models
for coastal aquifers (Eeman et al., 2012; Mollema and Antonellini, 2013; Oude Essink et al., 2010) and to integrate vadosezone processes more fully into groundwater models.
Vadose zone models transform weather data that are provided as input (precipitation and potential evapotranspiration) into
recharge fluxes arriving at the groundwater table. The choice ofthe temporal resolution of weather data (hourly, daily, monthly,
yearly) can significantly affect time distribution and values ofrecharge. This problem was investigated previously mainly
within the context of solute transport through the unsaturated
zone (e.g. Foussereau et al., 2001; Persson and Saifadeen, 2016;
Vero et al., 2014), with results showing that using longer
averaging periods (e.g., daily vs. hourly) of weather data leads
to less recharge and longer solute travel times. A systematic
study by Batalha et al. (2018) evaluated groundwater recharge
in different climatic regions of Brazil. They found that recharge
rates calculated using daily climatological data were up to 9
times larger than those obtained using yearly averages. On theother hand, saltwater intrusion studies are often performed
using monthly (e.g., Mollema and Antonellini, 2013; Prieto etal., 2006) or yearly (e.g., Oude Essink et al., 2010) weather
data. Thus, the question arises about the possible influence ofthe temporal resolution of weather data on modeling results.

Simulations of freshwater lens recharge and salt/freshwater interfaces using the HYDRUS and SWI2 packages for MODFLOW

It is generally accepted that comprehensive representationsof water flow in the vadose zone can be obtained using the
Richards equation (e.g., Healy, 2010), although this approach
generally does not account for some processes that may be
important at the local scale, such as preferential flow (e.g.,Šimůnek et al., 2003) or air entrapment (e.g., Szymańska et al.,2016). An extended form of the Richards equation can be used
to describe flow under both unsaturated and saturated conditions, thereby offering a unified three-dimension (3D) description of flow for the entire subsurface domain. Numerical codessuch as SUTRA (Voss and Provost, 2010), or HydroGeoSphere(Therrien et al., 2010) are based on this approach and include
the capability to simulate density-dependent flow essential forsaltwater intrusion modeling. However, due to high nonlinear
nature of the Richards equation and the need for fine spatial
discretizations close to the land surface, the requirements forcomputational power and time are often prohibitive for fully
coupled saturated-unsaturated 3D models.
Alternatively, less complex models can be used also to describe flow in the vadose zone, especially close to the land
surface. These models are often based on relatively simple
water balance principles (e.g., SWAT (Neitsch et al., 2011) or
HELP (Schroeder et al., 1994)) or kinematic wave routing (e.g.,the UZF1 package for MODFLOW (Niswonger et al., 2006)).
The usefulness of these type of models has been confirmed in a
number of studies (e.g., Chang et al., 2016; Holding and Allen,
2015; Luoma and Okkonen, 2014). They are typically very
efficient computationally. However, due to their simplified
nature, they may not be able to represent properly some
important features of unsaturated flow. For example, neither
UZF1 nor HELP allows for upward flow due to capillary action, while SWAT greatly simplifies water flow below the rootzone. Additionally, the kinematic wave model used in UZF1 isbased on the assumption of a homogeneous soil profile, which
is often not the case.
In addition to the above two groups of models, anotherpossible modeling approach of intermediate complexity would
be to couple a 3D (or 2D horizontal) model for flow in thesaturated zone with one or more 1D models for flow in representative unsaturated zone soil profiles, based on the Richards
equation. Such an intermediate approach would still allow for awide range of unsaturated zone processes, but with less computational overhead compared to the complete 3D variablysaturated flow models. The intermediate concept was
implemented in the HYDRUS package for MODFLOW (Beegum et al., 2018; Seo et al., 2007; Twarakavi et al., 2008),
which is based on the widely-used standalone numerical HYDRUS-1D code (Šimůnek et al., 2008ab, 2016). HYDRUS-1D
simulates the transient movement of water, solute, and heat invariably-saturated soil profiles using finite element discretization in space and fully implicit discretization in time. The model has been validated in a large number of applications(Šimůnek et al., 2008b, 2016), including several related torecharge (e.g., Scanlon et al., 2002; Šimůnek et al., 2016) andsaltwater intrusion (Werner and Lockington, 2004).
The HYDRUS package for MODFLOW was originally written for MODFLOW-2000 (Seo et al., 2007). Its evaluation andcomparison with other MODFLOW-related packages for vadose zone flow was presented by Twarakavi et al. (2008) andBailey et al. (2013). Applications to regional scale groundwater
flow, as well as comparisons with the UZF1 package, are described in Leterme et al. (2013, 2015). Currently, the package isbeing developed further in cooperation between researchgroups from the University of California Riverside, the Belgian
Nuclear Research Centre (SCK•CEN), and the Gdańsk Univer-sity of Technology (Beegum et al., 2018). The most recent
version of the HYDRUS package is compatible with MODFLOW-2005 (Harbaugh, 2005). So far, only a fraction of theprocesses implemented in the original standalone HYDRUS-1D
program has been ported to HYDRUS for MODFLOW. Much
potential exists for further enhancements, including features
such as contaminant transport, macropore flow, or improved
models of the hydraulic functions.
The SWI2 package (Bakker et al., 2013) is the latest releaseof the Seawater Intrusion (SWI) Package for MODFLOW. This
package allows simulations of variable-density groundwaterflow and seawater intrusion in multi-aquifer coastal systems.
Vertically integrated variable-density flow is based on the
Dupuit approximation in which an aquifer is vertically discretized into zones of differing densities, separated from eachother by defined surfaces representing interfaces (Bakker et al.,2013). The simplified, sharp-interface approach to densitydependent flow implemented in the SWI and SWI2 packages
allows for considerable savings of computer time, as comparedto solutions of the coupled variable-density flow and salt
transport equations. In fact, each aquifer can be represented byonly a single layer of cells. However, the sharp-interface model
neglects diffusion and dispersion processes, which may be
important in some circumstances. Numerical tests (Dausman etal., 2010) showed that the sharp-interface models may not berealistic when considerable dispersion occurs across the interface, when the anisotropy ratio is very small, or in some systems where inversion occurs and a significant amount of vertical fingering is present.
Our research reported in this contribution has two main
goals. The first objective is to carry out an evaluation of the
HYDRUS package for MODFLOW as a tool for simulating
recharge in conjunction with the SWI2 package (Bakker et al.,2013) for saltwater intrusion modeling. While HYDRUS-1D
has been used previously in studies related to saltwater intrusion, the HYDRUS package for MODFLOW has to the best of
our knowledge not yet been applied in this particular context.
The second objective is to investigate the influence of selectedfactors on the amount of recharge as well as evolutions of thegroundwater table and the saltwater-freshwater interface during
the simulations. These factors include uncertainty in soil hydraulic parameters (describing the retention and hydraulicconductivity functions of unsaturated soils), the presence ofvegetation, and the temporal resolution of weather data (precipitation and reference evapotranspiration).
MATERIALS AND METHODSHYDRUS code
Numerical simulations were performed using a modified
version of MODFLOW-OWHM (Hanson et al., 2014), which
includes both the SWI2 (Bakker et al., 2013) and HYDRUS
packages for MODFLOW (Twarakavi et al., 2008). We used
the updated version of HYDRUS for MODFLOW (Beegum etal., 2018), which represents an improvement over the original
version as described by Seo et al. (2007) and Twarakavi et al.
(2008), and which is compatible with MODFLOW-2005 and
MODFLOW-OWHM. Here, we provide a more detailed description of the current version of the HYDRUS package, while
details about MODFLOW-OWHM and SWI2 can be found in
Hanson et al. (2014) and Bakker et al. (2013), respectively.
The HYDRUS package performs simulations of verticalflow in the vadose zone overlying the MODFLOW model
domain. Simulations are carried out for a number of userdefined soil profiles. In principle, it is possible to associate a247

Adam Szymkiewicz et al.

separate soil profile with each cell of the MODFLOW horizontal grid, but this approach may lead to unacceptably long simulation times for larger problems. Thus, it is recommended todivide the flow domain into a number of regions in such a waythat each region encompasses several MODFLOW cells having
approximately the same vadose zone flow and transport parameters. The division into regions should be based on such criteria
as soil type, vegetation, depth to groundwater, and weather
data. It is also possible to define HYDRUS profiles only for apart of the MODFLOW cells for which a more detailed description of unsaturated flow is required, while using alternative
methods to specify recharge rates for the remaining cells.
The top of each HYDRUS profile corresponds with the soil
surface, while the bottom should be located within the uppermost aquifer defined in MODFLOW, below the expected watertable fluctuation zone. Each profile may consist of an arbitrary
number of soil layers having different hydraulic properties. The
number and position of nodes used to discretize each vadosezone profile are specified by the user. For reasons of numericalstability and efficiency, small nodal spacings (on the order of 1cm) close to the soil surface are recommended. Details of thenumerical scheme used to solve the Richards equation in HYDRUS are given by Šimůnek et al. (2008a).
Boundary conditions at the soil surface are implemented as
so-called atmospheric conditions, which can be described by
the following set of inequalities (Neuman et al., 1974; Šimůnek
et al., 2008a):
I max ≤ − K∂h− K ≤ E max∂zhdry ≤ h ≤ h pond

(1a)
(1b)

where Imax is the maximum possible infiltration rate (equal toprecipitation if interception is neglected), given as a negative
number, K is the soil hydraulic conductivity at the surface, h is
the pressure head, z is the vertical coordinate (oriented positively upwards), Emax is the maximum possible evaporation rate
(potential evaporation), hpond is the maximum ponding depth at
the soil surface, and hdry is the minimum allowed value of thepressure head at the surface. These values are provided by users
at an arbitrary temporal resolution, and not related to the stressperiods defined in MODFLOW. The type of the boundary
condition at the soil surface is automatically switched between
the flux and pressure head boundary conditions, depending
upon the prevailing weather and soil conditions. When thepressure head boundary condition is used, the potential flux is
reduced to the actual flux (e.g., actual infiltration and/or evaporation). Details about the numerical implementation of theatmospheric boundary conditions can be found in Šimůnek etal. (2008a).
Another important feature of the HYDRUS package is an
explicit representation of water uptake by plant roots. The rootzone depth and root density distribution is specified for eachsoil profile, according to available data on vegetation. The
maximum potential transpiration rates are specified by the user
at a temporal resolution matching the precipitation and evaporation data. At each HYDRUS time step, the current maximum
transpiration rate is distributed among the nodes in the root
zone, proportionally to the defined root density distribution, tocalculate the nodal potential root water uptake. Actual root
water uptake, appearing as a source term in the Richards equation, is obtained by multiplying potential root water uptake
using a stress response function, which in turn depends on the248

soil water pressure head. For our simulations we used a simplification of the well-known stress response function of Feddes etal. (1978), as given by Huang et al. (2011). Resulting values ofroot water uptake are either equal to or smaller than the potential root water uptake.
HYDRUS-MODFLOW coupling
The sequential coupling of MODFLOW and HYDRUS involves exchanging information about the recharge rates and theposition of the groundwater table. During each MODFLOW
time step, the HYDRUS package is called to solve the 1D
Richards equation for each defined soil profile. Due to thevariability and nonlinearity of flow in the unsaturated zone, the
HYDRUS solution is generally performed using much smaller
time steps than the MODFLOW solution. HYDRUS time steps
are adjusted automatically within a user-defined range according to the performance of the iterative solver. The initialcondition for each soil profile is specified by users for the first
MODFLOW time step. The atmospheric boundary condition at
the top of each profile is used to represent precipitation, irrigation and/or evaporation as described above. The boundary
condition at the bottom of each profile is given as the averagevalue of the groundwater head in the uppermost aquifer, calculated for all MODFLOW cells associated with a particular soil
profile.
The unsaturated flow simulation is performed for a time
period corresponding to the current MODFLOW time step.
Average flow through the bottom of the vadose zone profile
during a particular MODFLOW time step is then used as theaverage recharge (or capillary rise) rate for a particular time
step. This recharge rate is then included into the discretized
equations for the MODFLOW cells corresponding to a given
profile (in the same manner as recharge values defined in the
RCH package). Next, the solution for saturated zone flow is
obtained with MODFLOW for the new time level, thus providing new values of groundwater heads in the upper aquifer, withthe solution proceeding to the next time step. Pressure head andwater content distributions in each HYDRUS profile are updated next to reflect the new position of the water table and thevertical flux in the considered HYDRUS profile at the end of theprevious MODFLOW time step (Beegum et al., 2018).
Note that the HYDRUS package for MODFLOW does notconsider full coupling of saturated and unsaturated flow, suchas in the VSF (Thoms et al., 2006), VS2DI (Healy, 2008), or
SUTRA (Voss and Provost, 2010) codes. Rather, the HYDRUS
package is now merely a tool to facilitate and improve
estimations of infiltration and evapotranspiration fluxes, which
otherwise could be achieved using the standalone HYDRUS1D code or other codes solving the 1D Richards equation, as
implemented for example by Contractor and Jenson (2000) andKamps et al. (2008). The applied loose coupling scheme does
not allow complete consistency between the HYDRUS and
MODFLOW solutions due to the three reasons, First, during
each MODFLOW time step, the HYDRUS simulation is
performed with a constant groundwater head at the bottom of
the profile, which is subsequently updated instantaneously atthe beginning of a new time step. In reality, the groundwatertable varies continuously in time in response to vadose zone
flow and other processes such as pumping. Second, thegroundwater head at the bottom of each HYDRUS profile represents an average for all MODFLOW cells associated with aparticular soil profile, while in reality it is spatially variable.
And third, MODFLOW solves the groundwater flow equation
for an unconfined aquifer using a value of specific yield (Sy),Simulations of freshwater lens recharge and salt/freshwater interfaces using the HYDRUS and SWI2 packages for MODFLOW

which is defined by the user independently of the vadose zone
simulations. The value of Sy determines the change in the
groundwater table position in response to the specified recharge
flux. Consequently, the specific yield also determines the updated boundary condition at the bottom of each HYDRUS
profile during the next time step. We are aware that the value of
specific yield should be strongly related to the vadose zone
conditions above the water table; this issue is a subject of ongoing research.
The simplifications mentioned above result of the need tocompromise between accuracy and computational efficiency inthe proposed approach. As a consequence, it is not possible tocalculate a global mass balance for both the saturated and unsaturated zones. The mass balance should be considered separately for the MODFLOW groundwater domain (taking into
account the coupled fluxed) and for each HYDRUS profile ateach MODFLOW time step. In the simulations described below
the HYDRUS and MODFLOW mass balance errors for all time
steps were below 1%.
Simulation setup
Simulations were carried out for a 2D cross-section througha sandy peninsula along the Polish Baltic coast, as shown in
Fig. 1. Although the width of the peninsula was 600 m, thesolution domain extended offshore 200 m on each side. The
elevation of the ground surface varied from –5 m at the edgesof the domain to +5 m in the center of the peninsula, while thebottom of the aquifer was assumed to be horizontal at an elevation of –40 m. The MODFLOW grid consisted of a single rowof 100 cells in a single layer representing the unconfined aquifer. The vertical boundaries and the bottom of the domain wereassumed impermeable. A third-type boundary condition (ageneral head boundary, GHB) was prescribed on those parts ofthe soil surface that were below the sea level, similarly as in theexamples described by Bakker et al. (2013), using the GHB
package of MODFLOW. The total freshwater head, H = h+z,
was assumed to be 0 m for the GHB cells, while the conductivity was set to 100 m2/d. The SWI2 option ISOURCE was set to-2 for the GHB cells so that any water infiltrating into the aquifer was considered to be seawater, while water exfiltrating to
the sea had the salinity corresponding to the upper part of thecell. In all simulations, the number of salinity zones in the
SWI2 package was set to 2 (i.e., a single zeta-surface was defined). The density of freshwater (ρf) was 1 g/cm3, and thedensity of saltwater (ρs) was 1.0065 g/cm3, typical for south
Baltic sea water. The aquifer was assumed to be isotropic andhomogeneous with a saturated conductivity (Ks) of 10 m/day
and a specific yield (Sy) of 0.25.
The central (land) part of the aquifer was subject to recharge
from precipitation. Simulations of vadose zone flow were carried out using the HYDRUS package. Two regions with representative soil profiles were defined based on the average depth
to the groundwater (GW) table. Region 1 included two 200-m
wide strips of land parallel to the shoreline on each side of thepeninsula where the GW table was up to 2 m deep on average.
Region 2 consisted of the central 200 m of land where the GW
table was expected to be about 4 to 5 m deep. The soil profiles
were discretized using 41 nodes in Region 1 and 81 nodes inRegion 2. The vertical node spacing was non-uniform, ranging
from 1 cm close to the ground surface to 9 cm at the bottom of
the profile.
The simulations were performed for three types of sand,characterized by the van Genuchten-Mualem hydraulic functions of the following form:n −mθ ( h ) = θ R + (θ S − θ R ) 1 + (α h ) 


(2a)



n −1
n −m 

1 − (α h ) 1 + (α h )  K ( h) = KS KR ( h) = KS 
m/2
1 + (α h )n 



2

(2b)

where θ is the volumetric water content, θS is the water content
at field saturation, θR is the residual water content, h is the soil
water pressure head (negative in the unsaturated zone), K is thehydraulic conductivity, KS is the hydraulic conductivity at fieldsaturation, KR is the relative hydraulic conductivity, and α, nand m are fitting parameters (with m = 1–1/n). Values of thesoil hydraulic parameters are listed in Table 1. The first set,with a relatively high saturated conductivity (400 m/d), is denoted as medium sand in the database of VS2DTI and VS2DHI
codes (Hsieh et al., 1999) and was used e.g. by Dan et al.
(2012). While such a high hydraulic conductivity is not typicalfor vertical flow in the vadose zone (except macropore flow),
we include this set as an example of highly permeable material.
The second set, denoted here as fine sand, is reported by Carsel
and Parrish (1988) as an average set of parameters for sandaccording to the USDA textural division. Similarly, the third set
represents an average set of parameters for a loamy sand according to the USDA classification (Carsel and Parrish, 1988).
In the following section, we report results obtained for each ofthese three soil types. Note that the change in hydraulicFig. 1. Assumed flow domain and specified boundary conditions (CD = conductance).

249

Adam Szymkiewicz et al.

parameters affected only the HYDRUS simulations for thevadose zone. The values of specific yield and the horizontal
hydraulic conductivity in the aquifer were kept the same for all
simulations. Thus, we neglected the influence of unsaturated
flow on the specific yield, as explained above.
The initial condition for each simulation was obtained byrunning the MODFLOW/HYDRUS simulation for a period of
40 000 days with a constant value of the infiltration flux (equal
to the assumed recharge rate) imposed at the soil surface. The
infiltration flux was specified as 300 mm/yr, based on estimates
given by Dyck and Chardabellas (1963) (cited in Hölting and
Coldewey, 2013) for sands with a sparse plant cover (about
60% of yearly precipitation). A long warm-up period was usedbecause of the need to obtain the steady position of the interface with the SWI2 package, which generally requires morethan 100 years (see Examples 3 and 4 in the documentation of
SWI2). A warm-up period for the vadose zone simulations was
also necessary, but it is much shorter, on the order of a single
year. As a basic check, we compared the steady-state results
obtained with the HYDRUS package with a constant infiltration
rate to the results obtained with the RCH MODFLOW package.
This package allows one to input a user-specified rate of recharge at the water table. While the results are not shown here,
the agreement in terms of both the groundwater table elevationand the saltwater-freshwater interface position was very good(differences in the position of the groundwater table and the
fresh-salt interface in the center of the domain were less that 1.5
cm, while the elevations agreed with the Ghyben-Herzberg
formula (Verruijt, 1968)).
Simulations were performed for a 5-year period using the available weather data from Gdańsk (2011–2015), starting from theinitial steady-state conditions. Daily values of the reference
potential evapotranspiration (ET0) for the simulation period areshown in Figure 2, while the daily precipitation values areshown in Figs. 3, 5 and 7. Average annual precipitation during
the 5-year period was 550 mm/yr and the average annual ET0317 mm/yr. ET0 was estimated using the method of Grabarczyk, developed in Poland (Grabarczyk and Żarski, 1992),
which gives evapotranspiration values somewhat lower

Fig. 3. Daily values of recharge for fine sand without vegetation,
calculated using daily weather data; top: HYDRUS region 1, bottom: HYDRUS region 2.

Table 1. Parameters of the van Genuchten - Mualem hydraulicmodel for the soils used in this study.
Parameter

θRθS

KS [m/d]α [1/m]n = 3.08

Medium sand
0.022
0.375
400
4.31
3.10

Fine sand
0.045
0.430
10
14.50
2.68

Loamy sand
0.057
0.410
3.5
12.40
2.28

Fig. 4. Evolution of the water table position (top) and the freshwater/saltwater interface for different types of sand, using daily
weather data and assuming no vegetation.

Fig. 2. Daily reference evapotranspiration rates for the 5-year
simulation period.

250

than the widely used Penman-Monteith FAO equation. The
formula requires only daily values of the average air temperature and the average vapor pressure deficit.

Simulations of freshwater lens recharge and salt/freshwater interfaces using the HYDRUS and SWI2 packages for MODFLOW

Fig. 5. Daily values of recharge for fine sand with vegetation,
calculated using daily weather data in two HYDRUS regions; top:
HYDRUS region 1, bottom: HYDRUS region 2.

Fig. 6. Evolution of the water table position (top) and the freshwater/saltwater interface for different types of sand, using daily
weather data and assumed the presence of vegetation.

The total evapotranspiration flux includes both evaporationfrom the soil surface and transpiration of water taken up by
plant roots from deeper soil layers. Several studies (e.g., Comteat al., 2014; Kamps et al., 2008) reported that water uptake

Fig. 7. Daily values of recharge for fine sand without vegetation,
calculated using monthly weather data; top: HYDRUS region 1,
bottom: HYDRUS region 2.

by vegetation is an important factor influencing groundwater
resources of freshwater lenses. In order to account for this, we
performed two series of simulations. In the first series, we
neglected the presence of plants and assumed that ET0 was
equal to potential evaporation. In the second series, the reference evapotranspiration was distributed as a source term over
the root zone, i.e., ET0 corresponded to potential transpiration(interception was not considered). While both of these cases
were significant simplifications of reality, they can beconsidered as indicators for the possible range of variability ofrecharge for a varying soil cover by vegetation. We assumed
that the depth of the root zone was 2 m, with the root density
distribution decreasing linearly from a maximum value at thesoil surface to zero at a depth of 2 m. For the stress response
function, which describes the reduction in water uptake due todecreasing pressure heads in unsaturated soil, we used themodel of Huang et al. (2011), which is a simplified version of
the well-known formulation by Feddes et al. (1978). In thisapproach, root water uptake is assumed to be maximum forpressure heads larger than h1 (equal to -3 m in our study), zerofor pressure heads less than the wilting point (assumed at
h2 = –150 m), and changing linearly for pressure heads between
these two values. Such a model is considered representative ofconiferous forests (see Huang et al., 2011, and references therein). More details on the implementation of root water uptake in
HYDRUS-1D can be found in Šimůnek et al. (2008a). The two
scenarios represent typical conditions of coastal sandbars of theBaltic coast, which can be covered by bare dunes, grassland, orsparse pine forests.
The coupled codes were also used to investigate the effect of
the temporal resolution of weather data. We used either daily,
weekly, or monthly averages of the precipitation and evapotranspiration rates. Regardless of the temporal resolution of the251Adam Szymkiewicz et al.

weather data, we used a constant MODFLOW time step of 1
day and allowed the time steps in HYDRUS to vary within arange from 10–5 to 1 day. All simulations assumed a ponding
depth of zero and a minimum soil water pressure head at thesoil surface of –10000 m. The input data for the groundwatermodel were partly prepared using the ModelMuse interface(Winston, 2009). Since ModelMuse does not support HYDRUS
for MODFLOW, all data required by this package were prepared using spreadsheet software.
Sensitivity analysis
Additionally, a sensitivity analysis was performed bychanging values of the van Genuchten-Mualem parameters offine sand, one at a time, according to the following scheme: KS
was decreased or increased by one order of magnitude (KS = 1
or 100 m/d), α was changed by +/– 50% (α = 7.25 and 21.75
m–1), and n was changed by +/– 0.5 (n = 2.18 and 3.18). Such arange of parameter variability may reflect uncertainty in theunsaturated zone parameters as measured in the field or estimated from basic soil data via pedotransfer functions. For example, a study by Meyer et al. (1997) provides the following
range of variability for α and n in the USDA sand texturalgroup: KS = 0.3 m/d to 16 m/d, α = 6.87 to 22.6 m–1, and n =1.95 to 3.62. According to other authors, the hydraulic conductivity of sands may vary in a wider range. For example, Freeze
and Cherry (1979) gave values of KS from about 1 to 1000 m/d
for clean sands. These intervals are wider than those assumed inour study; combinations of hydraulic parameters that differ
from the above intervals may still be considered physicallyadmissible.

Table 2. Average annual recharge rates [mm/yr] obtained for
different simulation settings.
Scenario
Medium sand Fine sand Loamy sand
Daily weather data
476
449
416
No vegetation
Daily weather data
359
315
306
With vegetation
Weekly averaged weather data
370
354
341
No vegetation
Weekly averaged weather data
317
304
301
With vegetation
Monthly averaged weather data
338
326
318
No vegetation
Monthly averaged weather data
338
326
318
With vegetation

Table 3. Minimum (ζmin) and maximum (ζmax) elevations of the
freshwater-saltwater interface and the average annual recharge rate
for different values of the soil hydraulic parameters.
Scenario
Fine sand(Basic case)
Ks =100 m/d
Ks = 1 m/d
α = 29 /m
α = 7.25 /m
n = 3.08
n = 2.08

ζmin [m]ζmax [m]Average recharge[mm/yr]–28.96–28.45449–29.02–28.88–29.11–28.83–29.07–28.79–28.45–28.45–28.45–28.45–28.45–28.45

461
432
474
423
468
405

RESULTS AND DISCUSSION
Results obtained for the various scenarios described in theprevious section are presented in Tables 2 and 3 and Figures 3
to 8. The different scenarios are compared in terms of the average yearly recharge rate during the 5-year simulation period,
daily recharge values in the two HYDRUS regions, and evolutions of the groundwater table and the freshwater-saltwaterinterface in the center of the domain.
Influence of the soil hydraulic properties
The results of simulations performed using daily values ofprecipitation and reference evaporation for the three types ofsand are presented in Table 2 and Figures 3 and 4. For thesesimulations we neglected the presence of vegetation. Thismeans that the maximum evaporation flux at the soil surface
was set to ET0 as indicated in the previous section. Table 2
shows average annual values of recharge computed for the 5year period. They varied from 476 mm/yr for medium sand to
416 mm/yr for loamy sand, which corresponds to a relative
difference of about 13%. Figure 3 shows daily values of recharge in two regions, obtained using the fine sand hydraulicparameters. The distribution of recharge was found to be highlyvariable in time, generally following the pattern of precipitation. However, it can be seen that the depth of the vadose zone
has a visible influence on the temporal distribution of recharge
and its maximum daily values. In Region 2, where the depth togroundwater was larger than in Region 1, the recharge distribution was much smoother throughout the year, with peak values
largely reduced. On the other hand, recharge rates between precipitation periods were larger in Region 2 because of longer
redistribution time of soil moisture in deeper profiles.

252

Fig. 8. Evolution of the water table position (top) and the freshwater/saltwater interface for different types of sand, using monthly
weather data and assuming no vegetation.

As shown in Figure 4, the water table fluctuated following
periods of high and low recharge, with a maximum amplitude
of about 0.4 m. The freshwater-saltwater interface also showedseasonal fluctuations, with a general decreasing trend. The final
position of the interface was about 0.5 m lower than the initial
position. This can be explained by the fact that average yearlySimulations of freshwater lens recharge and salt/freshwater interfaces using the HYDRUS and SWI2 packages for MODFLOW

recharge rates during the simulation period (from 416 mm/yr to
476 mm/yr) were significantly higher that the recharge rateused obtain the initial steady-state condition (300 mm/yr). The
results shown in Figure 4 are similar for all three soil types,
except that for loamy sand the interface elevation decreased at alower rate than for medium and fine sand. We also note thatrecharge values were larger during the first 3 years of the simulations, resulting in a visible lower position of the interface. In
contrast, the last 2 years were characterized by less recharge
and more stable positions of the water table and the interface.
The simulated position of the interface was generally higher
than its elevation calculated using the Ghyben-Herzberg formula for a particular water table position. While the GhybenHerzberg approach assumes steady-state conditions, numericalsimulations represent transient conditions. A notable time lag
exists in the adjustment of the interface location for varying
positions of the water table.
Results of the sensitivity analysis are presented in Table 3.
The range of variability of average recharge rate and the final(minimum) position of the interface were consistent with therange of variability observed for the three sandy soils obtained
earlier. An increase in any parameter (i.e., KS, α, or n) produced
an increase in recharge and a decrease in the final position of
the interface. Modifying KS by an order of magnitude produceda change of less than 4% in the average yearly recharge. Changing the van Genuchten parameter α by a factor of two resulted
in variations in recharge of less than 6%. On the other hand,increasing or decreasing n by 19% produced a 4.2% increase
and a 9.9% decrease in the average recharge rate, respectively.
We hence conclude that, for our particular setting, the mostsensitive parameter was n since relatively small changes of its
value were found to have a significant effect on recharge.
Influence of water uptake by plant roots
The second series of simulations assumed that the evapotranspiration flux was distributed over a 2-m deep root zone.
Results of these simulations are summarized in Table 2. The
data show that root water uptake has a pronounced effect on
recharge, and consequently also on the saltwater-freshwaterinterface position. Average yearly recharge rates decreased tovalues between 306 and 351 mm/yr, corresponding to a reduction of 24 to 30% compared to scenarios without vegetation.
Relative differences in the average recharge rates of various
soils were slightly larger than those for corresponding cases
without vegetation (16% vs. 13%).
The reduction in recharge due to vegetation is also reflectedin the daily values of recharge shown in Figure 5. In contrast to
the previous case, the elevation of the freshwater-saltwaterinterface increased during the simulation period (Fig. 6) for all
three soils. The increase was particularly visible during the last
two years, characterized by smaller recharge rates, which isconsistent with the results for the scenarios without vegetation.
Also, in this case, one can see somewhat larger differences
between soil types in terms of the interface position. The highest final elevation of the interface was obtained for medium
sand.
Influence of the temporal resolution of weather data
In order to investigate the effects of the temporal resolutionof weather data, simulations were carried out using weekly andmonthly averaged precipitation and reference evapotranspiration rates. The resulting values of average annual recharge arereported in Table 2. For the bare soil scenarios with weekly andmonthly averaged meteorological data, the recharge rate decreased by about 20 and 27% , respectively. Daily recharge
values are shown in Figure 7 for the bare fine sandy soil andmonthly averaged weather data. The plots show that the recharge patterns are distinctly different from those in Figure 3.
As expected, fluctuations were smoothed to a large extent. The
episodic nature of recharge, reported in the literature (e.g., Huntet al., 2008; Smerdon et al., 2008), is lost when precipitation
values are averaged over longer time periods. Similar conclusions on the influence of time averaging of the weather data
were also presented by Batalha et al. (2018). These authors
found that using long-term averages of precipitation and evapotranspiration will lead to severe underestimates of the recharge
rates as compared to the use of daily data. Their results forsemi-arid region in Brazil differed by a factor of 9 between
simulations using daily and yearly averages of weather data. The
discrepancy was found to be particularly significant for coarsetextured soils, which corroborates the results of our study.
The influence of the temporal resolution of weather data wasfound to be significantly smaller for the scenarios with vegetation. Table 2 shows that the maximum differences in averageannual recharge were only about 13% for medium sand. In
contrast to the case of bare soil (no vegetation), using monthly
averaged weather data did lead to higher recharge rates thanweekly averages. While this could be explained by a more
complex pattern of water flow in the presence of a relatively
deep root zone, these interactions remain to be investigatedfurther. For the setting considered in this study, monthly averages of weather data produced nearly the same values of recharge for the cases with and without vegetation (differenceswere less than 1 mm/yr).
Table 2 also shows that the influence of soil hydraulic parameters on recharge decreases when weather data are averaged
over longer time periods. While relative differences in recharge
between the three types of sand for the daily data were about
13%, they were only about 8% and 6% for the weekly andmonthly averaged data, respectively. Sensitivity analysis withrespect to the hydraulic parameters KS, α, and n were also carried out for the weekly and monthly weather data. While details
are not further reported here, the obtained ranges of recharge
values were similar to the ranges shown in Table 2 for the three
types of sand. Similarly as for daily data, increasing KS, α, or n
produced a slightly higher recharge rates when weekly ormonthly weather data were used, with the n parameter seemingly being the most sensitive parameter.
CONCLUSIONS
Based on preliminary evaluations carried out in this study,
we consider the HYDRUS package for MODFLOW to be a
promising tool for integrating vadose zone processes in asaturated groundwater flow model accounting for saltwater
intrusion. This package allows calculations of spatially andtemporally variable recharge rates using the physically-based
Richards equation and an appropriate set of water retention andhydraulic conductivity functions characterizing vadose zone
flow processes, and accounting for interactions with groundwater.
Our simulations indicated that for the particular system considered in this study, changes in soil hydraulic parameters had amoderate influence on the final position of the groundwatertable and the freshwater-saltwater interface. This may explainedby the fact that we considered only sandy soils, which are characterized by low water holding capacities and relatively highhydraulic conductivities. The sensitivity analysis showed that253

Adam Szymkiewicz et al.

the recharge rates are affected more by the α and n parametersof the van Genuchten function than by the hydraulic conductivity, which for our study was not the limiting factor for water
infiltration. On the other hand, the results were influenced
more significantly by the temporal resolution of the weather
data and by the presence of vegetation (i.e., land use). The useof monthly and weekly averaged values of precipitation andpotential evaporation led to recharge rates that were more than
20% smaller than those obtained using daily values. Similarly,
when the potential evapotranspiration rate was distributed over
the root zone instead of being assigned to the soil surface (i.e.,vegetated versus unvegetated soil surfaces), the resulting recharge rates decreased by more than 24%. Differences in therecharge rate directly influenced the evolution of the saltwaterfreshwater interface. For the scenarios with daily weather dataand no vegetation, the position of the interface decreased withtime due to relatively large recharge (larger than the initialsteady-state value). The opposite trend was observed for scenarios with vegetation or for weekly and monthly averaged
precipitation and evaporation data, which produced reductionsin the recharge rate and upward displacement of the freshwatersaltwater interface.
Acknowledgements. This work has been supported by the National Science Centre, Poland, in the framework of project
2015/17/B/ST10/03233 "Groundwater recharge on outwash
plain". Weather data were obtained from a station operated by
the Department of Hydraulic Engineering, Faculty of Civil andEnvironmental Engineering, Gdańsk University of Technology.
The authors would like to thank two reviewers for their detailed and insightful remarks, which helped to improve themanuscript.
REFERENCES
Bailey, R.T., Morway, E.D., Niswonger, R.G., Gates, T.K.,
2013. Modeling variably saturated multispecies reactive
groundwater solute transport with MODFLOW-UZF and
RT3D. Ground Water, 51, 5, 752–761. DOI:
10.1002/ird.1699.
Bakker, M., Schaars, F., Hughes, J.D., Langevin, C.D.,
Dausman, A.M., 2013. Documentation of the seawater
intrusion (SWI2) package for MODFLOW. USGS
Numbered Series “Techniques and Methods” 6-A46. US
Geol. Survey, Reston, VA, USA.
Batalha, M.S., Barbosa, M.C., Faybishenko, B., van Genuchten,
M.Th., 2018. Effect of temporal averaging of meteorological
data on predictions of groundwater recharge. Journal of
Hydrology and Hydrodynamics, 66, 2, 143–152.
Beegum, S., Šimůnek, J., Szymkiewicz, A., Sudheer, K.P.,
Nambi, I.M., 2018. Implementation of solute transport in thevadose zone into the 'HYDRUS package for MODFLOW'.
Groundwater, (under review).
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability
distributions of soil water retention characteristics. Water
Resources Research, 24, 5, 755–769.
Chang, S.W., Nemec, K., Kalin, L., Clement, T.P., 2016.
Impacts of climate change and urbanization on groundwater
resources in a barrier island. Journal of Environmental
Engineering, D4016001.
Comte, J.C., Join, J.L., Banton, O., Nicolini, E., 2014.
Modelling the response of fresh groundwater to climate andvegetation changes in coral islands. Hydrogeology Journal,
22, 8, 1905–1920.

254

Contractor, D.N., Jenson, J.W., 2000. Simulated effect ofvadose infiltration on water levels in the Northern Guam
Lens Aquifer. Journal of Hydrology, 229, 3, 232–254.
Dan, H.C., Xin, P., Li, L., Li, L., Lockington, D., 2012.
Capillary effect on flow in the drainage layer of highway
pavement. Canadian Journal of Civil Engineering, 39, 6,
654–666.
Dausman, A., Langevin, C., Bakker M., Schaars, F., 2010. A
comparison between SWI and SEAWAT – the importanceof dispersion, inversion and vertical anisotropy. In:
Proceedings of the 21st Salt Water Intrusion Meeting,
Azores, Portugal.
De Louw, P.G., Eeman, S., Siemon, B., Voortman, B.R.,
Gunnink, J., Van Baaren, E.S., Oude Essink, G., 2011.
Shallow rainwater lenses in deltaic areas with saline
seepage. Hydrology and Earth System Sciences, 15, 3659–
3678.
Dyck, S., Chardabellas, P., 1963. Wege zur Ermittlung dernutzbaren Grundwasserreserven. Ber. Geol. Ges. DDR, 8,
245–262.
Eeman, S., Zee, S.V.D., Leijnse, A., De Louw, P.G.B., Maas,C., 2012. Response to recharge variation of thin rainwater
lenses and their mixing zone with underlying saline
groundwater. Hydrology and Earth System Sciences, 16, 10,3535–3549.
Eeman, S., De Louw, P.G.B., Van der Zee, S.E.A.T. M., 2017.
Cation exchange in a temporally fluctuating thin freshwaterlens on top of saline groundwater. Hydrogeology Journal,
25, 1, 223–241.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of
Field Water Use and Crop Yield. John Wiley & Sons, New
York, NY.
Foussereau, X., Graham, W.D., Akpoji, G.A., Destouni, G.,
Rao, P.S.C., 2001. Solute transport through a heterogeneous
coupled vadose‐saturated zone system with temporally
random rainfall. Water Resources Research, 37, 6, 1577–
1588.
Freeze, R.A., Cherry, J.A., 1979. Groundwater. Prentice Hall.
Grabarczyk, S., Żarski, J., 1992. Próba statystycznej
weryfikacji
niektórych
wzorówokreślających
ewapotranspirację potencjalną. [An attempt at statisticalverification of selected formulae estimating potential
evapotranspiration.]Zeszyty
NaukoweAkademiiTechniczno Rolniczej w Bydgoszczy (Rolnictwo). (inPolish.)Hanson, R., Boyce, S., Schmid, W., Hughes, J., Mehl, S.,
Leake, S., Maddock III, Th., Niswonger, R., 2014. OneWater Hydrologic Flow Model (MODFLOW–OWHM),
Techniques and Methods 6-A51. US Geological Survey,
Reston, VA, USA. Available from: http://dx. DOI:
org/10.3133/tm6A51.
Harbaugh, A.W., 2005. MODFLOW-2005, the US Geological
Survey modular ground-water model: the ground-water flow
process. US Department of the Interior, US Geological
Survey, Reston, VA, USA, pp. 6-A16.
Healy, R.W., 2008. Simulating water, solute, and heat transport
in the subsurface with the VS2DI software package. Vadose
Zone Journal, 7, 632–639.
Healy, R.W., 2010. Estimating Groundwater Recharge.
Cambridge University Press, Cambridge, UK, 245 p.
Holding, S., Allen, D.M., 2015. From days to decades:
numerical modelling of freshwater lens response to climatechange stressors on small low-lying islands. Hydrology and
Earth System Sciences, 19, 2, 933–949.

Simulations of freshwater lens recharge and salt/freshwater interfaces using the HYDRUS and SWI2 packages for MODFLOW

Hölting, B., Coldewey, W.G., 2013. Hydrogeologie:
Einführung in die allgemeine und angewandte
Hydrogeologie. Springer-Verlag.
Houben, G., Post, V.E.A., 2016. How long does the recovery ofa freshwater lens take after a massive saltwater inundation?
Experiences from the island of Baltrum, Germany, after the
1962 flood disaster. In: Proceedings of 24th Salt Water
Intrusion Meeting and the 4th Asia-Pacific Coastal AquiferManagement Meeting, 4 –8 July 2016, Cairns, Australia.
Hsieh, P.A., Wingle, A.W., Healy, R.W., 1999. VS2DI: A
graphical software package for simulating fluid flow andsolute or energy transport in variably saturated porous
media. USGS Water-Resources Investigation Report 994130. US Geological Survey, Reston, USA.
Huang, M., Barbour, S.L., Elshorbagy, A., Zettl, J.D., Si, B.C.,
2011. Water availability and forest growth in coarse-textured
soils. Canadian Journal of Soil Science, 91, 2, 199–210.
Hunt, R.J., Prudic, D.E., Walker, J.F., Anderson, M.P., 2008.
Importance of unsaturated zone flow for simulating recharge
in a humid climate. Ground Water, 46, 4, 551–560.
Illangasekare, T., Tyler, S.W., Clement, T.P, .Villholth, K.G.,
Perera, A.P.G.R.L., Obeysekera, J., Gunatilaka, A.,
Panabokke, C.R, Hyndman, D.W., Cunningham, K.J.,
Kaluarachchi, J.J., Yeh, W.W-G., van Genuchten, M.Th.,
Jensen, K., 2006. Impacts of the 2004 tsunami on
groundwater resources in Sri Lanka. Water Resour. Res., 42,
W05201. DOI: 10.1029/2006WR004876.
Jocson, J.M. U., Jenson, J.W., Contractor, D.N., 2002.
Recharge and aquifer response: northern Guam lens aquifer,Guam, Mariana Islands. Journal of Hydrology, 260, 1, 231–
254.
Kamps, P.W.J.T., Nienhuis, P., Witte, J.P.M., 2008. Effects ofclimate change on the water table in the coastal dunes of the
Amsterdam Water Supply. In: Proceedings MODFLOW.
Leterme, B., Gedeon, M., Jacques, D., 2013. Groundwater
recharge modeling of the Nete catchment (Belgium) using
the HYDRUS 1D - MODFLOW package. In: Šimůnek, J.,
M.Th. van Genuchten, Kodešová, R. (Eds.): Proc. of the 4thInternational Conference "HYDRUS Software Applications
to Subsurface Flow and Contaminant Transport Problems,
Prague, Czech Republic, pp. 235–244. ISBN: 978-80-2132380-3
Leterme, B., Gedeon, M., Laloy, E., Rogiers, B., 2015.
Unsaturated flow modeling with HYDRUS and UZF:
calibration and intercomparison. In: Proc. MODFLOW andMore 2015. Integrated GroundWater Modeling Center, May
31–June 3, 2015, Golden, CO.
Luoma, S., Okkonen, J., 2014. Impacts of future climate change
and Baltic Sea level rise on groundwater recharge,
groundwater levels, and surface leakage in the Hanko
aquifer in southern Finland. Water, 6, 12, 3671–3700.
Mahmoodzadeh, D., Ketabchi, H., Ataie-Ashtiani, B.,
Simmons, C.T., 2014. Conceptualization of a fresh
groundwater lens influenced by climate change: A modelingstudy of an arid-region island in the Persian Gulf, Iran.
Journal of Hydrology, 519, 399–413.
Meyer, P.D., Rockhold, M.L., Gee, G.W., 1997. Uncertainty
analyses of infiltration and subsurface flow and transport for
SDMP sites. Rep. NUREG/CR-6565, PNNL-11705. U.S.
Nuclear Regulatory Commission, Washington, DC.
Mollema, P.N., Antonellini, M., 2013. Seasonal variation innatural recharge of coastal aquifers. Hydrogeology Journal,
21, 4, 787–797.
Neitsch, S.L., Williams, J.R., Arnold, J.G., Kiniry, J.R., 2011.
Soil and water assessment tool theoretical documentation

version 2009. Texas Water Resources Institute, College
Station, TX.
Neuman, S.P., Feddes R.A., Bresler, E., 1974. Finite element
simulation of flow in saturated-unsaturated soils considering
water uptake by plants. Third Annual Report, Project No.
A10-SWC-77. Hydraulic Engineering Lab., Technion,
Haifa, Israel.
Niswonger, R.G., Prudic, D.E., Regan, R.S., 2006.
Documentation of the unsaturated-zone flow (UZF1)
package for modeling unsaturated flow between the land
surface and the water table with MODFLOW-2005 (No. 6A19).
Oude Essink, G.H.P., Van Baaren, E.S., De Louw, P.G., 2010.
Effects of climate change on coastal groundwater systems: amodeling study in the Netherlands. Water Resources
Research, 46, 10.
Persson, M., Saifadeen, A., 2016. Effects of hysteresis, rainfall
dynamics, and temporal resolution of rainfall input data insolute transport modelling in uncropped soil. Hydrological
Sciences Journal, 61, 5, 982–990.
Prieto, C., Kotronarou, A., Destouni, G., 2006. The influence oftemporal hydrological randomness on seawater intrusion incoastal aquifers. Journal of Hydrology, 330, 1, 285–300.
Sadurski, A., Borawska, J., Burczyk, T., 1987. Warunkihydrogeologiczne i hydrochemiczne Mierzei Helskiej[Hydrogeological and hydrochemical conditions of Hel
Peninsula.] Kwartalnik Geologiczny, 31, 4, 767–782. (InPolish.)Scanlon, B.R., Christman, M., Reedy, R.C., Porro, I., Šimůnek,J., Flerchinger, G.N., 2002. Intercode comparisons for
simulating water balance of surficial sediments in semiarid
regions. Water Resources Research, 38, 12, 1323. DOI:
10.1029/2001WR001233.
Schroeder, P.R., Dozier, T.S., Zappi, P.A., McEnroe, B.M.,
Sjostrom, J.W., Peton, R.L., 1994. The Hydrologic
Evaluation of Landfill Performance (HELP) Model:
Engineering Documentation for Version 3, EPA/600/R94/168b. US. Environmental Protection Agency, Risk
Reduction Engineering Laboratory, Cincinnati.
Seo, H.S., Šimůnek, J., Poeter, E., 2007. Documentation of the
HYDRUS Package for MODFLOW-2000, the U.S.
Geological Survey Modular Ground-Water Model, GWMI
2007-01. International Ground Water Modeling Center,
Colorado School of Mines, Golden, Colorado, 96 p.
Šimůnek, J., Jarvis, N. J., van Genuchten, M.Th., Gärdenäs, A.,
2003. Review and comparison of models for describing nonequilibrium and preferential flow and transport in the vadosezone. Journal of Hydrology, 272, 1, 14–35.
Šimůnek, J., Šejna, M., Saito, H., Sakai, M., van Genuchten,
M.Th., 2008a. The HYDRUS-1D software package forsimulating the one-dimensional movement of water, heat,and multiple solutes in variably-saturated media, Version
4.0. HYDRUS Software Series 3. Department ofEnvironmental Sciences, University of California, Riverside,
CA, USA, 315 p.
Šimůnek, J., van Genuchten, M.Th., Šejna M., 2008b.
Development and applications of the HYDRUS andSTANMOD software packages and related codes. Vadose
Zone Journal, 7, 2, 587–600.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2016. Recent
developments and applications of the HYDRUS computer
software packages. Vadose Zone Journal, 15, 7, 25 p. DOI:
10.2136/vzj2016.04.0033.
Sinclair, P., Galvis, S.C., Bosserelle, A.L., Post, V.E.A.,
Werner, A., 2016. Sustainability of freshwater lenses in atoll

255

Adam Szymkiewicz et al.

environments. In: Proceedings of 24th Salt Water Intrusion
Meeting and the 4th Asia-Pacific Coastal AquiferManagement Meeting, 4 –8 July 2016, Cairns, Australia.
Smerdon, B.D., Mendoza, C.A., Devito, K.J., 2008. Influenceof subhumid climate and water table depth on groundwater
recharge in shallow outwash aquifers. Water Resources
Research, 44, W08427.
Stuyfzand, P.J., 2016. Formation and hydrogeochemistry of a
freshwater lens on a sandbar island in saltwater lake
Grevelingen, Netherlands. In: Proceedings of 24th Salt Water
Intrusion Meeting and the 4th Asia-Pacific Coastal AquiferManagement Meeting, 4–8 July 2016, Cairns, Australia.
Sulzbacher, H., Wiederhold, H., Siemon, B., Grinat, M., Igel,
J., Burschil, T., Günther, T., Hinsby, K., 2012. Numericalmodelling of climate change impacts on freshwater lenses on
the North Sea Island of Borkum using hydrological andgeophysical methods. Hydrology and Earth System
Sciences, 16, 10, 3621–3643.
Szymańska, P., Tisler, W., Schütz, C., Szymkiewicz, A.,
Neuweiler, I., Helmig, R., 2016. Experimental andnumerical analysis of air trapping in a porous medium withcoarse textured inclusions. Acta Geophysica, 64, 6, 2487–
2509.
Therrien, R., McLaren, R. G., Sudicky, E. A., Panday, S. M.,
2010. HydroGeoSphere: A three-dimensional numericalmodel describing fully-integrated subsurface and surface
flow and solute transport. Groundwater Simulations Group,
University of Waterloo, Waterloo, ON.
Thoms, R.B., Johnson, R.L., Healy, R.W., 2006. User’s guideto the Variably Saturated Flow (VSF) process for
MODFLOW. Techniques and Methods 6-A18. US
Geological Survey, Reston, VA.
Trglavcnik, V., Robinson, C., Morrow, D., White, D., Paquin,
V., Weber, K., 2016. Effect of tides, waves and precipitation
on groundwater flow dynamics on Sable Island, Canada. In:
Proceedings of 24th Salt Water Intrusion Meeting and the 4th
Asia-Pacific Coastal Aquifer Management Meeting, 4–8
July 2016, Cairns, Australia.

256

Twarakavi, N.K.C., Šimůnek, J., Seo, H.S., 2008. Evaluatinginteractions between groundwater and vadose zone using
HYDRUS-based flow package for MODFLOW. Vadose
Zone Journal, 7, 2, 757–768.
Vandenbohede, A., Mollema, P.N., Greggio, N., Antonellini,M., 2014. Seasonal dynamic of a shallow freshwater lensdue to irrigation in the coastal plain of Ravenna, Italy.
Hydrogeology Journal, 22, 4, 893–909.
Vero, S.E., Ibrahim, T. G., Creamer, R. E., Grant, J., Healy, M.
G., Henry, T., Kramers, G., Richards, K.G., Fenton, O.
2014. Consequences of varied soil hydraulic andmeteorological complexity on unsaturated zone time lag
estimates. Journal of Contaminant Hydrology, 170, 53–67.
Verruijt, A., 1968. A note on the Ghyben-Herzbeg formula.
International Association of Scientific Hydrology Bulletin,
13, 4, 43–46. DOI: 10.1080/02626666809493624.
Voss, C.I., Provost, A.M., 2010. SUTRA: A model for
saturated-unsaturated, variable-density groundwater flow
with solute or energy transport. USGS Water- Resources
Investigations Report, 02-4231, U.S. Geological Survey,
Reston, VA
Werner, A.D., Lockinton, D.A., 2004. The potential for soil
salinization above aquifers influenced by seawater intrusion.
In: Proc. 13th International Soil Conservation Conference,
Brisbane, paper No. 790.
Werner, A.D., Bakker, M., Post, V.E., Vandenbohede, A., Lu,
C., Ataie-Ashtiani, B., Simmons, C.T., Barry, D.A., 2013.
Seawater
intrusion
processes,
investigation
andmanagement: recent advances and future challenges.
Advances in Water Resources, 51, 3–26.
Winston, R.B., 2009. ModelMuse: a graphical user interface for
MODFLOW-2005 and PHAST. US Geological Survey,
Reston, VA.
Received 19 June 2017
Accepted 10 January 2018