Journal of Hydrology 605 (2022) 127283

Contents lists available at ScienceDirect

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

Estimation of groundwater recharge in a shallow sandy aquifer using
unsaturated zone modeling and water table fluctuation method
Anna Gumuła-Kawęcka a, Beata Jaworska-Szulc a, Adam Szymkiewicz a, *,Wioletta Gorczewska-Langner a, Małgorzata Pruszkowska-Caceres a, Rafael Angulo-Jaramillo b,
Jirka Šimůnek cabc

Gdańsk University of Technology, Faculty of Civil and Environmental Engineering, Gdańsk, Poland
Univ Lyon, Université Claude Bernard Lyon 1, CNRS, ENTPE, UMR 5023 LEHNA, F-69518 Vaulx-en- Velin, France
University of California Riverside, CA, USA

A R T I C L E I N F O

A B S T R A C T

Keywords:
Groundwater recharge
Hydrus-1D
Water table fluctuation
Vadose zone flow

Quantification of groundwater recharge is one of the most important issues in hydrogeology, especially in view ofthe ongoing changes in climate and land use. In this study, we use numerical models of 1D vertical flow in thevadose zone and the water table fluctuation (WTF) analysis to investigate local-scale recharge of a shallow sandy
aquifer in the Brda outwash plain in northern Poland. We show that these two methods can be jointly used to
improve confidence in recharge estimation. A set of preliminary numerical simulations based on soil watercontent measurements from 4 grassland and pine forest profiles provided a wide range of recharge estimates(263 mm to 839 mm for a 3-year period). Additional simulations were performed with the lower boundary
condition specified as a functional relationship between the groundwater table elevation and the rate ofgroundwater outflow from the vertical profile (horizontal drains boundary condition). In this way, we could
reproduce the water table fluctuations resulting from recharge and lateral discharge to nearby lakes. The
agreement between simulated and observed groundwater levels differed depending on the specific set of pa­
rameters characterizing vadose zone flow, which allowed us to find the most representative parameter sets andrefine the range of plausible recharge estimates (501 mm to 573 mm per 3 years). The recharge rates from WTF(410 mm to 606 mm per 3 years) were in good agreement with numerical simulations, providing that the effectof the natural recession of groundwater table due to lateral outflow was considered (master recession curve
method). Our results show that: (i) the proposed approach combining 1D vadose zone modeling and WTF im­proves recharge estimation, (ii) multiple types of observations, including groundwater table positions, areneeded to calibrate and validate vadose zone flow models, and (iii) extended periods of observations and sim­
ulations are necessary to capture year-to-year variability in the recharge rates.

1. Introduction
Sustainable management of groundwater resources requires an ac­
curate estimation of groundwater recharge flux. This knowledge isessential to quantify the safe yield of an aquifer in the hydrologicalbalance of catchment (Bouwer, 1989; de Vries & Simmers, 2002;
Sophocleous, 1991). The groundwater recharge flux also plays a sig­
nificant role in the transport of contaminants because it determines thetravel time through the unsaturated zone (e.g., Sousa et al., 2013;
Szymkiewicz et al., 2018b; 2019), contaminant load reaching theaquifer (e.g. Szymkiewicz et al., 2020), directions of groundwater flowin shallow aquifers (e.g., Scanlon et al., 2002a), and consequently
vulnerability of aquifers, and in particular groundwater intakes, tohealth-related hazards. This is particularly true for shallow, unconfined
aquifers, such as aquifers located on outwash plains. An outwash plain,also named ’sandur’, is mostly a sandy deposit formed by glacial abla­tion water. It is a geomorphological form widely occurring in the glaciallandscape of the North European Plain, Fennoscandia, northern Siberia,
Canada, and the northern US. A significant spatial extent and high hy­
draulic conductivity of the outwash deposits make sandurs abundant
groundwater bodies commonly used to supply water.
Several studies have been carried out in recent years to quantify* Corresponding author.
E-mail address: adams@pg.edu.pl (A. Szymkiewicz).
https://doi.org/10.1016/j.jhydrol.2021.127283
Received 1 October 2021; Received in revised form 24 November 2021; Accepted 28 November 2021
Available online 11 December 2021
0022-1694/© 2021 Published by Elsevier B.V.

A. Gumuła-Kawęcka et al.

Journal of Hydrology 605 (2022) 127283

groundwater recharge on outwash plains and other sandy deposits ofsimilar character in temperate and cold climates. Their results are
summarized in Table 1 in terms of the average annual recharge and therecharge/precipitation ratio ω. In colder areas, such as Finland, Sweden,
and western Canada, where melting of snow and ice rapidly reaches
groundwater (Şen, 2015), recharge can reach almost 70% of the yearly
precipitation. In warmer and more humid regions, recharge is usually
lower because of less snow and longer vegetation period, which increase
evapotranspiration (ET). Groundwater recharge in the North European
Plain (e.g., Poland or western Russia), with a more continental climate(precipitation 500 – 700 mm/yr), is about 9 – 24% of precipitation inforested or agricultural areas. In the lowlands with more abundant
rainfall (greater than700 mm/yr), such as northern Germany, Belgium,or the eastern coast of North America, recharge ranges between 13 and
34% of precipitation in agricultural areas, 24 – 51% in coniferous for­ests, and 28 – 53% in grasslands.
A variety of methods is used to estimate groundwater recharge. They
can be broadly categorized as physical techniques, water-budget
methods, numerical modeling, and tracer methods (Healy, 2010; Scan­lon et al., 2002a). Within these groups, we can distinguish methods
based on measurements from different hydrological compartments, i.e.,
the unsaturated zone, the saturated zone, and surface water. Using
several methods is recommended to obtain reliable results (Healy, 2010;
Lerner et al., 1990; Scanlon et al., 2002a; Simmers, 1997). In the studies
of Callahan et al. (2012), Huet et al. (2016), Jie et al. (2011), Krogulec(2010), Liu et al. (2014), Szilagyi et al. (2011), and Zomlot et al. (2015),
recharge estimations in temperate climate were conducted with at least
two methods.
Numerical models of vadose zone flow are increasingly applied toobtain recharge estimates, including most of the studies reported inTable 1. Unsaturated zone models can be based on the solution of the
Richards equation (e.g., HYDRUS-1D (Šimůnek et al., 2013)), the kine­
matic wave equation (e.g., MODFLOW UZF (Niswonger et al., 2006)), or
soil water balance equations (e.g., WETSPASS (Batelaan & De Smedt,
2001) and HELP (Schroeder et al., 1994)). The required computer pro­grams are often freely available. They can represent a number of factors
influencing recharge in a physically based manner, which makes theman ideal tool for capturing the variability of recharge in space and timeand for investigating scenarios of land use or climate change. However,such models require a significant number of parameters that may bedifficult to obtain. Consequently, the resulting recharge estimates areprone to uncertainties resulting from (a) spatial variability of soil (e.g.,layering), (b) uncertainty in the water retention and hydraulic con­ductivity functions, and (c) uncertainty in parameters determining root
water uptake and ET (e.g., leaf area index LAI and root depth). In orderto provide a meaningful prediction of recharge, vadose zone models
must be calibrated and validated against field measurements and othermethods of recharge estimation (Jie et al., 2011). This is not a simple
task since data on vadose zone flow and ET are difficult to obtain.
Moreover, including this data in the calibration process does not guar­antee that a unique set of model parameters is derived. Brunner et al.
(2012) showed that at least in some scenarios using ET data in cali­
bration constrains the vadose zone model better than soil moisturemeasurements, while the most valuable input for model calibration isrepresented by groundwater heads. Jie et al. (2011) proposed an
approach to cross-validate recharge rates from the soil water balance
model using observations of groundwater table fluctuations. They
calculated specific yield based on the recharge estimate from the soilmodel and the increment of groundwater table observed in the corre­
sponding time period. The obtained values of specific yield wereTable 1
Groundwater recharge based on selected studies of unconfined sandy aquifers.
Reference

Location

Land cover

Applied method

Precipitation[mm/yr]Recharge[mm/yr]Recharge/
Precipitation ratioω[-]
Cherkauer & Ansari(2005)
Nastev et al. (2005)
Smerdon et al. (2008)

Southeastern Wisconsin,
USASouthwestern Quebec,
CanadaNorthern Alberta, Canada

Wendland et al. (2008)

Lower Saxony, Germany

Mixed crops and pasture

Dripps & Bradbury(2010)
Krogulec (2010)

Northern Wisconsin, USA

Grass
Coniferous forest
Forest-agricultural

Grinevskii &Novoselova (2011)
Leterme & Mallants(2011)
Leterme et al., (2012)

Moscow Artesian Basin,
western Russia
Dessel, Belgium

Forest
Field
Coniferous forest

Dessel, Belgium

Grass

Kurylyk & MacQuarrie(2013)
Ala-aho et al. (2014)

New Brunswick, eastern
CanadaNorthern Finland

Coniferous anddeciduous forest
Pine forest

Graf & Przybyłek(2014)

Poznań Plateau, western
PolandForestWu (2014)

Tärnsjö, central Sweden

Coniferous forests

Pozdniakov et al.
(2015)
Åberg et al. (2019)

Voronezh, southwestern
RussiaNorthern Finland

Forest1
2

Vistula River valley, central
PolandAgricultural, urbanizedBaseflow analysis

989

224–248

0.23–0.25

Agricultural, urbanized

Water table fluctuations

930–1130

131–225

0.13–0.221

Boreal forest

Numerical modeling (inhouse code)
Numerical modeling(GROWA)
Soil – water balance model

410

45

0.11

740

150–2502

0.20–0.34

570–980

Groundwater model(MODFLOW)
Water table fluctuations
Unsaturated flow model(HYDRUS-1D)
Unsaturated flow model(HYDRUS-1D)
Unsaturated flow model(HYDRUS 1D)
Soil-water balance model(HELP3)
Numerical modeling(CoupModel)
Soil-water balance model(WetSpass)589

160–435
140–425
56–95

0.28–0.53
0.24–0.51
0.10–0.16

657
667
899

75
138
100
239

0.13
0.21
0.15
0.27

899

314

0.35

1230

541(average)

0.44 (average)591

363 (average)

0.61(average)

517–579
541(average)600

46–13181(average)

0.09–0.24
0.15 (average)

210–3100.35–0.52

558–672

130

0.19–0.23560

315–378

0.56–0.68

ForestUnsaturated flow model(CoupModel)
Unsaturated flow model(HYDRUS-1D)
Groundwater model(MODFLOW)

compared to the average precipitation of 1030 mm/yr.
locally more than 250 mm/yr.
2

A. Gumuła-Kawęcka et al.

Journal of Hydrology 605 (2022) 127283

compared to the expected values for the aquifer to check if the measured
water table fluctuations are consistent with the simulated recharge
episodes.
Groundwater level measurements are simpler to perform and usuallymore available than vadose zone data or ET measurements. However,including them in a state-of-the-art 1D vertical model of vadose zone
flow based on Richards’ equation is not straightforward, as it is related
to the choice of the boundary condition (BC) at the bottom of the soilcolumn. In vadose zone flow simulations, the bottom BC is often
formulated in terms of the water pressure head corresponding to themeasured position of the water table, which can be constant or variablein time. If the water table elevation varies rapidly in time, using it as the
BC can lead to unphysical water fluxes at the bottom of the unsaturated
zone, which causes errors in recharge estimations (Beegum et al., 2018).
Moreover, using a pressure head-based BC does not allow reproducing
the natural sequence of water table fluctuations resulting from infiltra­tion and ET processes occurring in the vadose zone. Alternatively, thebottom of the soil column can be assumed impermeable (Brunner et al.,2012). In this case, one can represent water table fluctuations due tovadose zone processes, but the influence of lateral groundwater flow isneglected, which is unrealistic for many aquifers, especially those dis­charging to nearby surface water bodies. The third option is the freedrainage boundary condition, where the bottom of the soil column isabove the groundwater table, and outflow is driven by gravity, with thevertical pressure head gradient equal to zero (e.g., He et al., 2021; Healy,
2010; Ines & Mohanty, 2008; Šimůnek et al., 2013; Zendehboudi et al.,2012). The bottom flux gives an estimate of recharge, but the exact
arrival time of recharge at the water table cannot be reproduced, andinformation about water table elevation is not included in the model.
In this study, we suggest using yet another option for the bottom BC
in vadose zone simulations, i.e., the pressure head-dependent flux. The
bottom of the soil profile is below the groundwater table, and theoutflow rate from the saturated zone is calculated as a function of thecurrent water table position. The function contains two parameters,
which are calibrated using the observed water table levels. The cali­
bration procedure can be performed for different sets of vadose zone andET parameters to find which parameters allow for the best agreement
between the observed and calculated water table fluctuations. In thisway, uncertainty related to the estimation of vadose zone parameters
and the resulting recharge rates is reduced.
Besides their utility for the calibration and validation of vadose zone
flow models, groundwater level observations in shallow unconfined
aquifers can be used to estimate recharge by the water table fluctuation
method (WTF) (Healy, 2010; Scanlon et al., 2002a). This method re­quires only one parameter – the specific yield Sy. Unfortunately, Sy isdifficult to determine, as it depends on the type of soil, the depth to thegroundwater table, and water dynamics in the vadose zone (Healy,2010). The simplest variant of WTF is the RISE method (Rutledge, un­published manuscript, after Delin et al., 2007), which assumes thatrecharge events only occur when there is an increase in the groundwaterlevel. The RISE method tends to underestimate recharge. In reality, thechanges in the groundwater level are driven by other processes besides
recharge and ET, especially lateral flow towards streams, lakes, andother surface water bodies. Thus, recharge also often occurs when thewater table is stagnant or decreases, but this is masked by lateral
outflow. In order to account for this factor, the MRC method has been
proposed (Heppner & Nimmo, 2005, Delin et al., 2007; Heppner et al.,2007, Jie et al., 2011). In this approach, a master recession curve (MRC)
is identified for a specific observation well using observations from dry
periods. MRC represents a decrease in the groundwater level due tolateral outflow and other discharge processes in the absence of recharge.
The rate of the water table decline (ΔH/Δt) decreases with a decreasingwater table elevation H, and several approaches have been proposed toestablish a functional relationship between these two quantities (Nimmo
et al., 2015). The water level increments are calculated with respect tothe MRC, leading to higher recharge estimates than using the RISE

method (Delin et al., 2007). Thus, the results of the WTF analysis areaffected by uncertainty related to the choice between the RISE and MRC
methods, the estimation of MRC, and specific yield.
The aim of our study was to combine vadose zone flow modelingbased on the Richards equation with pressure head-dependent flux as
the bottom BC and the WTF analysis, in order to reduce the uncertainty
in recharge estimation. The two methods are applied to local-scale
recharge estimates on the experimental site on the Brda outwash plain
in northern Poland. Another objective was to assess time variability of
the recharge in a 3-year observation period, which contained an
extremely wet year, followed by an extremely dry year. The influence ofland cover (grassland vs. pine forest) and the related estimates of LAI
and root depth on groundwater recharge were also investigated.
2. Study area
Field measurements were carried out on a site near Tuchola town innorthern Poland. This area is a part of the Brda sandur, considered to be
the largest outwash plain in Poland. The experimental site is almost
completely surrounded by four lakes (Fig. 1). The land is covered mainly
by pine forests, with some open grasslands, crop fields, and dispersed

Fig. 1. Location of the observation site and soil profiles (modified from
Gumuła-Kawęcka et al., 2018). H, L and A denote the highest, the lowest andthe average elevation of the water table in the 3-year observation period.
3

A. Gumuła-Kawęcka et al.

Journal of Hydrology 605 (2022) 127283

low buildings (summer houses and farms). According to the data fromthe weather station in Chojnice (35 km from the experimental site) forthe period 1951–2019, the yearly average air temperature in the regionis 7.4 ◦ C, and the average annual precipitation is 581 mm.
Figure 1 presents the soil profiles obtained from 4 observation
boreholes (P1 and P4 in pine forest, P2 and P3 on grassland). In eachprofile, we can distinguish two sand layers (Sa1 and Sa2), separated by aloamy sand / sandy loam layer (SaLo), occurring about 1.5 m below theground level and 1.5 to 2 m thick. Moreover, in P1 and P4, we distin­
guished a layer of topsoil (TS) based on the observed differences in soil
water dynamics with respect to Sa1. Below the vadose zone is a shallow
aquifer, with the water table at a depth of 6 to 7 m and the bottom at a
depth of about 15 m. The aquifer is composed of fine/medium sand withgravel. It is recharged by infiltrating precipitation and drained by sur­rounding lakes which are connected to Brda river with a stream. Below
the shallow aquifer is a layer of weakly permeable glacial till (15 to 24 m
below ground level (b.g.l.), not shown in the profiles), which overlies adeeper, confined aquifer. The lower aquifer, also composed of glacio­
fluvial sand with gravel, is recharged by the upper aquifer, which was
confirmed by Oficjalska and Gregosiewicz (2000) and Muter (2002).
Measurements performed in another borehole in October 2016 showedthat the piezometric level in the deeper aquifer was about 0.6 m below
the groundwater table in the shallow aquifer, which corresponds to avertical hydraulic gradient of 0.095.
There is a regional tendency of lowering groundwater tables inshallow unconfined aquifers within the Brda outwash plain (JaworskaSzulc et al., 2017). The analysis of groundwater levels in 10 monitoring
wells for the period of 2003–2016, obtained from the monitoring
network, showed a long-term trend of water table lowering since 1990.
In most cases, the groundwater table has been below the monthly
average water level since 2006. The tendency is related to lower infil­
tration recharge, likely caused by climate change (less snow cover,irregular and more intensive rainfall events in summer). The lowering of
the water table in the shallow aquifer may reduce or eliminate therecharge gradient between the shallow and deeper aquifers, which could
be a long-term and potentially irreversible effect. Preliminary numericalsimulations of vadose zone flow based on literature data (GumułaKawęcka et al., 2018) showed significant year-to-year variability ofgroundwater recharge, following changes in annual precipitation in the
period 2003–2016. Regional-scale estimations of average annual
groundwater recharge in the Brda river watershed (from baseflow
measurements) were given by Hobot et al. (2012). They range from 155
mm/yr for the upper part of the watershed (above Tuchola) to 193 mm/
yr for the whole watershed, which corresponds to 26% and 33% of theaverage precipitation in Chojnice, respectively.

3.2. Soil hydraulic characteristics
Several methods were applied to estimate hydraulic characteristics
for soils at the site. They included pedotransfer functions based onparticle soil distribution, single ring infiltrometer tests, tension infil­
trometer tests, saturated hydraulic conductivity measurements with the
Aardvark permeameter (SoilMoisture Equipment Corp.), and measure­ments of soil water retention functions with sand-kaolinite apparatus. A
detailed presentation of the results is outside the scope of this report andwill be the subject of a separate study. Here we use data derived fromfield measurements to constrain calibration of vadose zone flow models.
We assumed that each soil layer is homogeneous and isotropic withineach specific profile, but parameters of the same layer may differ be­
tween profiles. Hydraulic characteristics of soils were described with thevan Genuchten-Mualem model (van Genuchten, 1980):
θ − θr= [1 + (α|h| )ng ]− mg
θs − θr(1a)[() mg ] 2
K = Ks Kr = Ks Seτ 1 − 1 − Se1/mg

(1b)

Se =

where Se is the effective saturation [-], θr is the residual volumetric watercontent, θs is the saturated volumetric water content, α is the parameter
related to the average pore size [L− 1], ng, mg are parameters related tothe pore size distribution [-], mg = 1–1/ng, K is the unsaturated hydraulicconductivity [LT− 1], Ks is the saturated hydraulic conductivity [LT− 1],
Kr is the relative hydraulic conductivity [-] and τ is a parameter related
to pore geometry [-].
The saturated water contents θs for the TS, Sa1, and SaLo layers wereset equal to the average porosity, as measured on undisturbed soil
samples. Since no undisturbed samples were taken from Sa2, we
assumed that it is more compact due to the self-weight of the soil and set
θs = 0.35, which was the minimum value measured in the upper layer
Sa1. The residual water contents θr were set to be slightly below thesmallest values measured by soil moisture probes in the TS, Sa1, and Sa2
layers, while in SaLo, we set it to 0.057 (an average value for loamy
sand, according to Carsel and Parrish (1988)). The saturated hydraulicconductivities Ks for the TS and Sa1 layers were obtained from single
ring infiltrometer tests and measurements with the Aardvark per­meameter and from permeameter measurements and pumping tests for
the Sa2 layer.
Pumping tests were conducted in piezometers to determine thehorizontal hydraulic conductivity of the aquifer. Wells were pumped
with a constant rate of about 4.6 m3/h. The water level was measured in
15-minute intervals until drawdown stabilized. Hydraulic conductivity
was calculated based on the Dupuit method with the Forchheimer
correction for partial penetration. The radius of drawdown was esti­
mated using the Kusakin empirical formula (Pazdro & Kozerski, 1990).
Due to these simplifying assumptions, the calculated hydraulic con­ductivity must be considered as approximate. Nevertheless, it is in good
agreement with the measurements using the Aardvark permeameter inthe same sand layer Sa2, about 1.5–2 m above the water table. The poreconnectivity parameter τ was set to the default value of 0.5 in all soil
layers, except Sa2 in P1, where τ was calibrated, to ensure better
agreement between field measurements of water content and simula­tions. The remaining parameters, i.e., α and n for all layers and Ks for
SaLo, were found by calibration of the numerical model described
below. Permeameter measurements of Ks for the SaLo layer were notincluded in the calibration due to the large variability of results andpotential errors.

3. Methods
3.1. Field measurements of hydrological data
Continuous soil and groundwater data measurements were carriedout from April 15, 2017, to April 14, 2020, in 4 observation profiles (P1P4). The water table level and the electric conductivity of groundwater
were recorded using sensors (baro-divers Decagon Devices CTD-10)
placed in piezometers, screened in the interval of approximately 5 to
11 m below the ground level (b.g.l.) in each piezometer. Volumetricwater content was measured in each profile using automatic sensors(TDR, Decagon Devices, GS3 and FDR, Delta_T Devices SM150) at four
depths (ranging from 20 cm to 450 cm b.g.l., Fig. 1).
Meteorological data (precipitation, temperature, air pressure, rela­tive humidity, wind speed and direction, and solar radiation) wascollected by a weather station located in an open field (P2) from April
15, 2017, to July 30, 2019. An additional pluviometer was placed in theforest (P1) to estimate interception. All measurements of weather, soil,
and groundwater data were recorded by automatic data loggers in 10minute intervals.

3.3. Numerical simulations with a constant water table
Numerical models of vadose zone flow were developed for eachprofile using the HYDRUS-1D computer program (Šimůnek et al., 2013).
HYDRUS-1D solves the Richards equation describing vertical flow in a4

A. Gumuła-Kawęcka et al.

variably-saturated soil profile.
()∂θ(h) ∂
∂h
∂k(h)
=
K(h)
+
− S(h)∂z∂t∂z∂z

Journal of Hydrology 605 (2022) 127283

were relatively mild (see section 4.1 and Table 2). For the same reason,
we did not consider soil freezing (preliminary simulations showed thatthe difference in recharge estimates with and without soil freezing is
negligible). According to de Beer’s law, PET was divided into the po­tential evaporation flux at the soil surface and potential transpiration in
the root zone (Šimůnek et al., 2013). We assumed nonlinear root dis­tribution with depth, according to Jackson et al. (1996):(2)

where θ is the volumetric water content [L3L− 3], t is time [T], h is thewater pressure head (negative in the unsaturated zone), z is the spatial
coordinate [L], K(h) is the hydraulic conductivity function of the un­saturated medium [LT− 1], and S(h) is a sink function representing water
uptake by plant roots [L3L− 3T− 1].
For spatial discretization of Eq. (4), HYDRUS-1D uses Galerkin-typelinear finite elements with mass lumping, while a fully implicit firstorder scheme is used for time discretization. Further details about thenumerical scheme can be found in the code documentation (Šimůnek
et al., 2013). The code has been validated in a large number of appli­
cations (Šimůnek et al., 2008; 2016), many of them related to the esti­
mation of groundwater recharge (e.g., Batalha et al., 2018; Grinevskii &Novoselova, 2011; Leterme et al., 2012; Leterme & Mallants, 2011;
Pozdniakov et al., 2015; Scanlon et al., 2002b; Šimůnek et al., 2016;
Szymkiewicz et al., 2018a; Twarakavi et al., 2008).
Two series of simulations were carried out. In the first series, theprofile depth was set to the average depth to groundwater tablemeasured in the observation period (as shown in Fig. 1), and the bottom
boundary condition was given as h = 0 = const. These simulations wereused to calibrate the models and estimate groundwater recharge. We
used uniform spatial discretization with 1 cm node spacing. In order toavoid the influence of the initial condition on the results, we applied a 3year warm-up period, starting with the water pressure distribution close
to field capacity and using weather data from the Chojnice station.
An atmospheric boundary condition was specified on the soil surface
using precipitation and potential evapotranspiration fluxes. No water
ponding was assumed to occur at the soil surface (instantaneous runoff),
and the minimum allowed pressure head on the soil surface in dry pe­
riods (separating stage 1 from stage 2 of evaporation) was set to − 1000
m. We used precipitation measured by rain gauge at P2 (grassland) inthe main simulation period from April 15, 2017, to April 14, 2020, andmeasurements from the Chojnice weather station in the warm-up period(April 15, 2014 –April 14, 2017). Interception was included according to
the model implemented in HYDRUS-1D, based on von Hoyningen-Hüne(1983), Braden (1985), and van Dam et al. (1997).
⎛⎞⎜
I = a∙LAI⎝1 −1⎟⎠
1 + P∙SCFa∙LAI

Y = 1 − βd(4)where Y is the cumulative root fraction (from 0 to 1), β = 0.943 forgrassland (P2 and P3) or 0.976 for temperate coniferous forest (P1 and
P4), and d is a depth in cm. The resulting root density sharply decreases
with depth, which is consistent with the fact that even in pine forests, themajority of roots typically occur within the uppermost 20 cm of the soil,
and only very few roots can be found below 1.5 m (Rutkowski et al.
2017; Sainju & Good, 1993). Root water uptake of grass and pine forest
was estimated with the Feddes et al. (1978) macroscopic model. The
plant-specific stress response function parameters were taken from De
Silva et al. (2008).
For each observation profile, we simulated two scenarios using thelower and upper estimates of LAI. In each scenario, the soil parameters
not derived from field measurements were found using inverse modelingby fitting the numerical solution to measured water contents in eachprofile. Time series of daily-averaged water contents in 30-days intervals
from all four depths (P1, P2, P4 – 20, 50, 150, 450 cm, and P3 – 10, 90,
180, 280 cm) were used for that purpose. The optimization of soil hy­
draulic parameters was carried out using the Marquardt-Levenberg al­
gorithm implemented in HYDRUS-1D. The recharge estimate
corresponding to the calibrated soil parameters was taken as the final
result of the simulation in each scenario.
3.4. Numerical simulations with lateral groundwater outflow
The second series of simulations aimed to validate the numericalmodels described in the previous section by reproducing observed var­iations in groundwater table elevations. For each profile, we used the setof calibrated parameters from the earlier simulations. The only changes
introduced in the second series of simulations were: (i) extending eachprofile to a depth of 8 m (keeping the same node spacing of 1 cm), to
include the upper part of the aquifer in the Sa2 layer, and (ii) setting thebottom boundary condition to a relationship between the hydraulichead and flux, to represent lateral groundwater discharge to nearby
lakes. For this purpose, we used the “Horizontal drains” feature of
Hydrus-1D. The simplest analytical formula describing flow to drain
available in Hydrus is a simplification of Hooghoudt’s equation:(3)

where I is interception [LT− 1], a is the interception constant [LT− 1], LAI
is leaf area index [-], P is precipitation [LT− 1], SCF is the soil cover
fraction [-], estimated as SCF = 1-exp(-0.463∙LAI). The interception
constant a was set to 0.25 mm for grassland (a typical value for crops)
and 0.37 mm for the forest. The latter value ensured the best agreement
between the interception calculated in the HYDRUS model and theinterception measured on-site (difference in precipitation recorded bygauges on the open field and in the forest).
Potential evapotranspiration (PET) was calculated with the PenmanMonteith equation (Allen et al., 1998), as implemented in Hydrus-1D,using daily weather data measured at the site (2017–20) and in Choj­
nice (2014–17). When calculating PET, we used a typical range of LAI
based on Scurlock et al. (2001). The average LAI in temperate evergreen
needle leaf forests is 5.47, with a standard deviation of 3.37. Since pine
forests in Poland tend to have sparse canopies (Jagodziński & Kałucka,
2008), we took LAI = 2.10 and 5.47 as our site’s lower and upper esti­mates, respectively. For grasslands, Scurlock et al. (2001) reports
average LAI = 1.71 (+/-1.19). As the area has a relatively dense grasscover, we took LAI = 1.71 and 2.90 as the upper and lower estimates,
respectively. LAI was assumed constant throughout the year since bothtypes of plant cover are perennial, and winters in the simulated period

qdr =

4Kdr
(H − Hdr )2 = C(H − Hdr )2
L2dr

(5)

where qdr is the lateral outflow (discharge to drains) per unit area, Kdr isthe horizontal hydraulic conductivity, Ldr is the drain spacing, H is thewater table elevation in the soil profile, Hdr is the water table elevationin the drain, and C=(4Kdr/Ldr2) is a lumped parameter representing
hydraulic conductance of the aquifer and aquifer/lake interface. Un­
fortunately, no detailed water level measurements in the lakes were
performed, and it likely changed during the observation period. On theother hand, groundwater flow was influenced by factors not represented
in Eq. (5), e.g., vertical hydraulic gradients or low permeability of thelake bottom sediments. In principle, it is possible to use more complex
models instead of Eq. (5), but then more parameters need to be esti­
mated or calibrated. For the sake of simplicity, we used Eq. (5) andconsidered Hdr and C as fitting parameters. Their values were estimatedby manual calibration, based on the difference between the observedand calculated positions of the water table in each profile, calculated as
the sum of square errors SSQ = Σ (Hobs - Hsim)2, where Hobs and Hsim are
the observed and simulated elevations of the groundwater table.
5

A. Gumuła-Kawęcka et al.

Journal of Hydrology 605 (2022) 127283

Table 2
Climate data for the experimental site: (1) from Chojnice weather station, (2) measured on site.
Period15 April 2014–14 April
2015
15 April 2015–14 April
2016
15 April 2016–14 April
2017
15 April 2017–14 April
2018
15 April 2018–14 April
2019
15 April 2019–14 April
2020

Average daily
temperature [◦ C]

Potential
evapotranspiration[mm/yr]Annual precipitation[mm/yr]Number of days with averagetemperature < 0 ◦ C

Number of days withsnow cover

LAI =
1.71

LAI =
5.74

9.3(1)

681.4

1034.6507.7(1)

42(1)

24(1)

9.1(1)

725.0

1106.8418.8(1)

36(1)

20(1)

8.8(1)

642.2

964.7

703.6(1)

51(1)

52(1)

8.2(2)

521.3

610.1

919.6(2)

54(2)

32(1)

10.2(2)

633.4

722.1

463.0(2)

42(2)

12(1)

10.1(2)

591.4

676.5

670.4(2)

9(2)

4(1)

HYDRUS-1D allows two options to incorporate the calculated flux qdr
in the discrete form of the Richards equation. It can be added to thebottom node as a boundary condition in the strict sense, or it can be
distributed over all nodes below the groundwater table as a sink term.
The latter option seemed more consistent physically with the assump­tion of horizontal groundwater flow and was used in our simulations.

above the highest position of the water table in each profile.
3.6. Water table fluctuation (WTF) analysis
The WTF analysis was based on water level measurements in pie­zometers, averaged for each day of the observation period. The daily
amount of recharge Ri is calculated as:

3.5. Estimation of specific yield

Ri = Sy ΔHi

We combined several approaches to derive a plausible range of
specific yield values for our site, including literature data (Johnson,
1967, Healy, 2010) and the empirical formula of Bieciński (Pazdro &Kozerski, 1990) commonly used in Poland, which links Sy to the hy­
draulic conductivity of the aquifer Ks:√̅̅̅̅̅
Sy = 0.117 7 Ks
(6)where Sy is the specific yield and ΔHi (greater than 0) is an increment in
the groundwater table elevation caused by recharge occurring between
day i-1 and i. A range of Sy values was obtained using the methods
described in the previous section. We applied two methods (RISE and
MRC) to calculate the increments ΔHi. The RISE method assumes thatrecharge only occurs if there is an increase in groundwater table
elevation between two subsequent days, so ΔHi = Hi – Hi-1 if Hi > Hi-1,
otherwise ΔHi = 0.
The MRC method assumes that in the absence of recharge, the watertable decreases daily by a specific amount ΔHMRCi, following the socalled master recession curve (MRC), which represents in a simplified
manner the discharge processes in the aquifer, in particular the lateral
outflow to nearby surface water bodies.
The groundwater table tends to decrease at a rate depending on its
current position, i.e., the higher the water table, the faster the lowering
rate (Heppner & Nimmo, 2005). MRC provides a functional relationshipbetween a daily decrement of the water table ΔHMRCi and the water table
elevation Hi-1 in periods without recharge. In this study, we used thesimplest, linear MRC:

where Ks is given in [m/d]. In the above formula, we used the estimateof Ks in the Sa2 layer in each profile. Finally, we also estimated Sy fromHYDRUS-1D simulations with a lateral outflow boundary condition. For
each day of the simulation, we calculated the volume of water in thelower section of the profile, which contained fluctuating water table.
The zone included in the calculations stretched from the bottom of the
profile to some height above the water table, as discussed later.
Assuming that the soil profile represents a column of 1 m2 cross-section,
we express the volume of water as the equivalent water height Hw inmeters. The specific yield Syi for a given day i was calculated as:
Syi =

ΔHwi
ΔHi

(8)

(7)

ΔHMRCi = A∙Hi− 1 + B

where ΔHwi is the daily change in water volume and ΔHi is the dailychange in groundwater table position (both quantities can be either
positive for increasing water table or negative for decreasing watertable). The final value of Sy was obtained as an average of Syi for all daysof the simulation period. Due to transient effects, one can obtain Syi
negative or larger than soil porosity for some days. These values were
excluded from the calculation of the average Sy.
Some consideration is needed to choose the upper limit of the zone
affected by the water table fluctuations, which is the basis for calculating
Hw. If the upper limit is too close to the water table, it is in the capillaryzone, and the resulting Sy is very small. If it is too close to the soil sur­face, Sy is affected by surface infiltration and evapotranspiration, which
are unrelated to water table fluctuations. A preliminary analysis forprofile P1 showed that the results were relatively consistent if the upper
limit of the calculation zone was set in the range from 1 m above thehighest position of the water table in the profile to the bottom of the
SaLo layer. Thus, we set the upper limit of the zone for calculating Sy 1 m(9)The coefficients A and B were fitted separately for each piezometer.
The fit was based on data from periods of continuous groundwater level
decrease longer than two weeks. More refined methods to estimate MRC
are available (Heppner & Nimmo, 2005; Jie et al., 2011), but weconsidered the approach described here as sufficiently accurate for the
purposes of our study. The daily increment of the water level due torecharge was then calculated as: ΔHi = Hi - Hi-1 + ΔHMRCi if Hi > (Hi-1ΔHMRCi), otherwise ΔHi = 0 (the decrement ΔHMRCi taken as a positive
number).
4. Results and discussion
4.1. Weather data and measured groundwater levels
Figure 2 shows daily precipitation totals measured at the P2 obser­
vation point from April 15, 2017, to April 15, 2020. There is a significant
6

A. Gumuła-Kawęcka et al.

Journal of Hydrology 605 (2022) 127283

Fig. 2. Precipitation and potential evapotranspiration at the observation site.

Fig. 3. Measured and simulated water contents in the forest profiles P1 and P4 (sensor depth: 20 cm (black), 50 cm (dark blue), 150 cm (green), 450 cm (light blue)).
7

A. Gumuła-Kawęcka et al.

Journal of Hydrology 605 (2022) 127283

difference between meteorological conditions during the three years ofobservation (Table 2). The one-year period starting on April 15, 2017,
was exceptionally wet – the amount of precipitation measured by a rain
gauge at profile P2 (open field) reached 920 mm (more than 150% of thelong-term average), with a mean temperature equal to 8.2 ◦ C. In
contrast, the one-year period starting on April 15, 2018, was unusuallywarm and dry, with an average temperature of 10.2 ◦ C and precipitationof 462 mm. The third year (April 15, 2019 – April 14, 2020) was warm
(10.1 ◦ C), but the rainfall total was slightly above average – 670.4 mm/
y. Among the 36 months of the observation period, 32 months had
monthly temperatures higher than the long-term average for theconsidered month, based on long-term (1951–2020) observation datafrom Chojnice station. Ten months in the observation period were
extremely warm, within 5% of the warmest months in the years
1951–2020.
Yearly precipitation measured by the rain gauge in the forest (P1)
was 687 mm in the first year, 321 mm in the second year, and 501 mm in
the third year, indicating that average interception by the forest canopy
was about 27% of yearly precipitation. This is compatible with the range
observed in pine woods in Poland (Puchalski & Prusinkiewicz, 1990).
Maximum and minimum daily values of PET calculated with thePenman-Monteith equation are also presented in Fig. 2. The annual PET
was the lowest in the first year (starting on April 15, 2017) – between
521.3 and 610.1 mm, while PET reached the highest values in the second
year (April 15, 2018 –April 14, 2019) – from 633.4 mm to 722.1 mm. In
the third year (starting on April 15, 2019), PET ranged between 591.4mm and 676.5 mm. The differences are caused mainly by different
weather conditions during vegetative periods in 2017, 2018, and 2019.
In 2017, the vegetation period (May – October) was relatively cold andrainy, with an average temperature of 14.5 ◦ C and a precipitation
amount of 682 mm. In contrast, in 2018, the vegetative period was
extraordinarily warm and dry, with an average temperature of about
16.5 ◦ C and a rainfall total of 234 mm. More mild conditions wereduring the vegetation period in 2019 – the mean temperature reached15.5 ◦ C, but the precipitation sum was 438 mm.
The measured values of volumetric soil water content are presented
in Figs. 3 and 4. A distinct pattern of seasonal variability can be seen,
with wet periods: Apr 2017 - Apr 2018, Oct 2018-Apr 2019, Oct 2019Apr 2020, and dry periods in Apr-Oct 2018 and Apr-Oct 2019. In thefirst year of observations, there is an exceptionally long wet period withhigh water contents at all sensors, which is not repeated in the subse­
quent drier years. Under the sandy loam layer, water content

Fig. 4. Measured and simulated water contents in the grassland profiles P2 and P3 (sensor depth: 2P –20 cm (black), 50 cm (dark blue), 150 cm (green), 450 cm(light blue), and 3P – 10 cm (black), 90 cm (dark blue), 180 cm (green), 280 cm (light blue)).
8

A. Gumuła-Kawęcka et al.

Journal of Hydrology 605 (2022) 127283

fluctuations in the Sa2 layer are generally much smaller than in the TS
and Sa1 layers. In the SaLo layer, the water content is significantly
higher than in other layers due to finer soil texture. The fluctuations inthis layer are larger in forest profiles than in grassland profiles, probably
due to water uptake by deeper tree roots in dry periods.
A clear variability pattern can also be seen in water table elevations
measured in four profiles (Fig. 5). In the first (wet) year, the ground­water table increased monotonically. Then it decreased, reaching thelowest elevation in September 2019, and rose again in the last months ofthe observations. During the long period of decline, we can distinguishan interval with a relatively low decline rate in autumn/winter 2018/19, corresponding to increased soil water contents in Figs. 3 and 4. In
contrast, during the vegetation periods in 2018 and 2019 (Apr-Oct), we
observe low soil water contents and a sharp groundwater table decline.
There are some breaks in water table measurements in P1 and P2 due tomalfunctions of the sensors. Also, the water table elevation in P2 isabove the other three piezometers starting from mid-2018 which isprobably caused by technical issues with the barodiver and the breaks inoperation mentioned above. Different flow dynamics might have also
played a role since P2 is the well farthest away from the lakes. Overall,
the soil moisture and groundwater table show strong seasonal and yearto-year variability related to large differences in annual precipitationand ET.

in the SaLo layer (obtained from inverse numerical solution) shows
significant differences between profiles, with the lowest values in P1 andthe highest in P2 (different by a factor of about 500). In most cases, thecalibrated values of α and n were sensitive to the choice of LAI, whichstrongly influenced evapotranspiration. However, we obtained a similar
fit quality in each of the eight scenarios, reported as RMSE in Table 2. In
terms of RMSE, the best fit was obtained for P2, and this is confirmed bya visual comparison of simulated and measured water contents in Figs. 3
and 4. On the other hand, in P2, the pattern of water content evolution isdifferent from other profiles since the water content is almost constant in
the SaLo layer, while in other locations, it decreased significantly during
the dry periods.
4.3. Recharge estimates from numerical simulations
Figure 6 shows annual and total (3-year) recharge values obtainedusing the simulations with a constant water table position and the final
calibrated sets of parameters in each scenario. They differ significantly
between the profiles and strongly depend on LAI, ranging from 263 mm
in P4 with an LAI of 5.47 to 839 mm in P2 with an LAI of 1.71 in thethree years. For the lower LAI estimates, recharge can be even more than
two times larger than for the upper LAI estimate due to smaller inter­ception and ET. There is a reasonable similarity in predictions for thetwo forest profiles (P1 and P4). There is no overlap in the results for thetwo grass profiles, with the lower estimate in P2 minimally exceeding
the upper estimate in P3. In contrast, the range calculated for P3 partly
overlaps with the estimates for P1 and P4.
There are also clear differences in the time distribution of therecharge flux arriving at the groundwater table. Four examples areshown in Fig. 7. In all cases, recharge occurs mostly during the initial(wet) year of the observations. However, in P1-LAI2.10 and P2-LAI2.90,
the flux increases and later decreases gradually over longer time periods,

4.2. Soil hydraulic characteristics and soil water contents
Parameters of the van Genuchten model for all profiles and soil layers
are summarized in Table 3. The saturated hydraulic conductivity ob­tained from field measurements in the Sa1 and Sa2 layers does not differ
substantially between the profiles, with Ks in Sa2 consistently higher
than in Sa1. In P1 and P4, the topsoil is characterized by a lower hy­
draulic conductivity than the underlying Sa1 sand layer. In contrast, Ks

Fig. 5. Groundwater table elevation measured in piezometers at four observation sites.
9A. Gumuła-Kawęcka et al.

Journal of Hydrology 605 (2022) 127283

by either Hdr or Ldr.
The agreement between the best fit and observations varied
depending on the soil profile and LAI. The lowest SSQ values (1.5⋅104)
were obtained for P1 with an LAI of 2.10 and P3 with an LAI of 2.90. Thisis confirmed by a visual inspection of the plots, which shows a relatively
good agreement between observations and the model for most of thesimulation period. For P4 with an LAI of 2.10, the SSQ is only slightly
larger (1.6⋅104), but the fit is clearly worse for the rising branch of thehydrogram, corresponding to the initial wet period. The largest SSQ
were obtained for P2, where the groundwater table measurements were
probably inaccurate, as discussed earlier, while on the other hand, therecharge estimates were higher than for other profiles. For the forest
profiles (P1 and P4), the fit was better for the lower LAI estimate (2.10),
while for the grass profiles (P2 and P3), it was better for the higher LAI
estimate (2.90). The results presented in Figs. 8 and 9 indicate thatsimulations for P1 and P4 with an LAI of 2.10 and P3 with an LAI of 2.90
provided recharge estimates, which are most representative for our sites,
i.e., they are in the best agreement with the observed pattern ofgroundwater level fluctuations. Consequently, total recharge in theobserved period can be estimated to be between 501 and 573 mm, asignificantly narrower range than the one resulting from the first seriesof simulations. The corresponding LAI falls between 2.10 and 2.90,
representing the upper estimates for grassland and lower estimates forpine forests.
The fitted values of conductance coefficient C fall in the range from1.32⋅10− 5 to 4.44⋅10− 5 cm.d− 1. The corresponding theoretical distance
between the observed profile and the lake (half of the drain spacing) can
be calculated from Eq. (5), assuming the aquifer permeability Kdr =3716 cm/d (an average value for the Sa2 layer from Table 2). The ob­tained distance is between 91 and 168 m. This is in reasonable agree­
ment with the real distances: 30 m for P1, 170 m for P2, 70 m for P3, and
90 m for P4. Note, however, that there is no clear relationship between
the fitted C value and the real distance to the lake at a specific point.

Table 3
Parameters of the soil hydraulic functions (values obtained using inverse
modeling are denoted by superscripts: 1 – simulations with a lower LAI estimate,2 – simulations with an upper LAI estimate).
Topsoil (TS)Sand (Sa1)

Sandy loam
(SaLo)

Sand Sa2

RMSE [-]

θr [-]
θs [-]
Ks[cm/d]
α [1/
cm]
n [-]
τ [-]
θr [-]
θs [-]
Ks[cm/d]
α[1/
cm]
n [-]
τ [-]
θr [-]
θs [-]
Ks[cm/d]
α [1/
cm]
n [-]
τ [-]
θr [-]
θs [-]
Ks[cm/d]
α [1/
cm]
n [-]
τ [-]

P1

P2

P3

P4

0.025
0.400
238
0.0071,2
2.761,
2.462
0.5

n/a

n/a

0.025
0.400
1356
0.0101,2
2.481,
2.322
0.5

0.045
0.400
813
0.0371,0.01921.701,
1.742
0.5
0.057
0.410
1021, 1072
0.0041,0.00621.801,
1.392
0.5
0.045
0.350
3767
0.0381,0.02923.451,
3.372
0.5

0.045
0.400
506
0.0141,2
1.631,
1.622
0.5

0.025
0.400
91
0.0181,0.01222.201,
1.882
0.5
0.025
0.400
502
0.0531,0.0172
1.291,2
0.5

0.057
0.410
351, 332
0.0091,0.0082
1.381,
1.392
0.5
0.025
0.350
3145
0.0371,0.0402
1.891,
1.872
0.5

0.057
0.410
0.21, 220.0051,0.00224.001,
1.432
0.5
0.025
0.350
4216
0.0361,
0.02122.991,
3.642
0.5

0.0141,0.01520.0271,0.02620.0231,
0.02220.057
0.410
0.21,20.0021,0.0012
3.611,
2.892
0.5
0.025
0.350
3737
0.0151,0.0142
3.531,
3.202
2.081,
1.972
0.0201,
0.02224.5. Estimation of specific yield

while in P2 and P4, separate recharge events are clearly marked as peak
flux values. Moreover, in P2, relatively large fluxes occur in the second
and third years of the observation period. In P1 and P3, there is only
small recharge in the second year and somewhat larger recharge in thelast year. Finally, in P4, there is practically no recharge in the second
year and only a small flux in the last year. These results are consistent
with the differences in soil hydraulic parameters and LAI between thesimulations.
Considering the whole 3-year period of observations, the recharge/
precipitation ratio varies from 13% to 41%, depending on the consid­
ered scenario. Note that despite such a wide variability, the R/P ratio
falls within the range of values reported in Table 1. Since we obtained areasonably good fit to the measured water content values in all cases,
additional information was required to reduce the uncertainty resulting
from the variability of soil characteristics and the lack of precise data on
evapotranspiration. Thus, another series of simulations was carried out
to reproduce the groundwater table’s observed fluctuations and narrow
the range of recharge predictions.

Healy (2010), following Johnson (1967), reports Sy = 0.21 for finesands with the standard deviation +/- 0.07. Since Ks estimates for Sa2
are quite similar in each profile (Table 3), the empirical formula ofBieciński, Eq. (6), gives a narrow range of Sy estimates: 0.191 to 0.200.
Sy derived from the second series of numerical simulations (with afluctuating water table) ranged from 0.229 (P3, LAI = 1.71) to 0.300(P4, LAI = 2.10), with an average of 0.269. These results are consistent
with the upper range of estimates given by Healy (2010), which can be
explained by the high uniformity of grain sizes in Sa2. Ultimately, in the
WTF analysis, we used 0.200 and 0.280 as the lower and upper estimatesof Sy for our sites.
4.6. Recharge estimations using WTF
The master recession curves for each profile are shown in Fig. 11.
There is a considerable scatter between the points representing therelationship between the groundwater table elevation (H) and a daily
decrement of the groundwater level (ΔHMRC). In each case, it waspossible to fit a linear trend, showing an expected decrease of the ab­solute value of ΔHMRC with decreasing H, but the results must be
considered as highly approximate.
Figure 12 presents groundwater table increments calculated using
the RISE and MRC methods for P3. Both methods produce a similar
pattern of recharge variability. However, in the MRC approach, there are
significantly more days with recharge, and the estimated recharge rates
are higher than in the RISE method. According to the MRC procedure,
recharge also occurs when there is a decline in the groundwater table,
but it is smaller than the decline resulting from the master recession
curve. This can be seen especially in the period from Oct 2018 to Apr
2019, characterized by a small decrease in the groundwater table, as

4.4. Water table fluctuations from numerical simulations
The water table fluctuations obtained using numerical simulations
with a lower boundary condition representing lateral outflow are shownin Figs. 8 and 9. In each case, the best fit resulting from manual cali­
bration of boundary condition parameters is presented. The results were
found to be quite sensitive to both Hdr and C, which facilitated theircalibration. Hdr determined the average elevation of the water table,
while C influenced both the average elevation and amplitude of fluc­tuations (Fig. 10). However, the shape of the curve is determined by the
distribution of recharge over time, and it is not fundamentally affected
10

A. Gumuła-Kawęcka et al.

Journal of Hydrology 605 (2022) 127283

Fig. 6. Recharge estimates obtained using different methods. For numerical simulations, the range of results corresponds to the assumed range of LAI values in eachprofile. For WTF, the range of results corresponds to the assumed range of specific yield values in the aquifer.

discharge to lakes is partly compensated for by recharge frominfiltration.
Recharge rates calculated for each profile using the RISE and MRC
methods and different estimates of Sy are shown in Fig. 6, together withthe results of numerical simulations discussed earlier. Several conclu­
sions can be drawn from this comparison. The estimates obtained using
the MRC method are more than twice as large as those obtained using
the RISE method, emphasizing the importance of including a recession
curve in the analysis. In most cases, the results of numerical simulations
are closer to MRC than to RISE, except for the second (driest) observa­
tion year, where the opposite is true. In this period, recharge is low andthus sensitive to inaccuracies in the approximation of the recession
curve.
Considering the whole 3-year period, we can notice a good agree­
ment between the MRC method and the HYDRUS-1D simulations for theforest profiles (P1 and P4) with the lower LAI (2.10) and for the P3
profile (grass with an LAI of 2.90). The numerical results fall between
the MRC estimates for low and high Sy (410 mm to 606 mm). In contrast,
the largest discrepancy between numerical modeling and the WTF
analysis occurs for P2. In this case, even for the highest value of specificyield, the MRC estimate is significantly smaller than the lowest value
obtained using the numerical model. These observations further support
the conclusions from Section 4.4 that the simulations for P1 with a low
LAI and P3 with a high LAI yield recharge patterns, which are the mostconsistent with the observed groundwater level fluctuations, while thelargest incompatibility is observed for P2. The disagreement for P2 isprobably caused by two overlapping effects: increased local recharge
due to the high permeability of the sandy loam layer and a possible

malfunction of the water level sensor.
Further cross-validation of the WTF method and numerical modelingcan be done following Jie et al. (2011). The sum of the water table in­crements calculated using the MRC approach for P1 is 2080 mm in the 3year period, while the corresponding recharge obtained using the
HYDRUS-1D simulation with an LAI of 2.10 is 571 mm. This provides anestimate of Sy (=571 mm/2080 mm) of 0.274, very close to the value
calculated using the water balance in the numerical solution with lateral
outflow (0.278) and within the range given by Healy (2010). Similarly,
for P3, the sum of increments is 2050 mm, and the recharge estimate
using Hydrus (LAI = 2.90) is 501 mm for Sy = 0.244. Again, this value isvery close to the estimate derived directly from the numerical simulation
with lateral outflow (0.233). This analysis increases the confidence inthe recharge estimates provided by numerical simulations.
4.7. Time variability of recharge and a recharge/precipitation ratio
The presented results show that groundwater recharge is a timevariable process and, at least on a short time scale, cannot be consid­
ered as steady-state flow. The time distribution of groundwater recharge
is strongly determined by the intensity and distribution of precipitation.
Abundant rainfalls in 2017 resulted in high inflow to the water table inautumn/winter 2017/2018. On the contrary, the drought in 2018
caused a drop in the water table and low percolation from the vadosezone to the aquifer in winter/spring 2019. In the last year of the ob­servations (2019/20), both precipitation and recharge were between thetwo earlier extreme years. In Table 4, we report annual and total
recharge estimates based on three numerical simulations that were11

A. Gumuła-Kawęcka et al.

Journal of Hydrology 605 (2022) 127283

Fig. 7. The recharge flux (vBot in Hydrus 1D) at the groundwater table obtained using numerical simulations with a constant water table position.

found to be the most representative for the study area (P1-LAI 2.10, P3LAI 2.90, P4-LAI 2.10). The ratio of annual recharge between the first(wet) and second (dry) year of observations ranged from 3.1 (P4-LAI
2.10) to 7.5 (P3-LAI 2.90). The ratios obtained using WTF were smaller,
but still significant: 2.3 to 3.6 in the RISE method and 1.8 to 2.7 in the
MRC method.
The time of the water table response to a single rain event depends on
the thickness of the unsaturated zone, soil type, and vegetation. Since
the magnitude of the water table rise is determined by the water content
in the unsaturated zone, the amount of precipitation during the pre­ceding weeks plays a significant role (Dripps, 2012). That is probably thereason why groundwater recharge following the storm at the end of June
2017, which occurred after wet months, is so apparent. On July 17,
2018, even more abundant rainfalls (47 mm/day) took place. However,this event did not result in any response of the water table, probably
because of the dry preceding spring and summer.
Some studies provide guidelines to estimate annual or multi-year
recharge in Central Europe as a specific fraction of precipitation,
depending on the soil type, land cover, and other factors (Dyck andChardebellas, 1963, cited after Hölting and Coldewey, 2019; Duda et al.,2011). Hölting and Coldewey (2019) provide an empirical formula forthe annual recharge-precipitation relationship based on an earlier study

of Dyck and Chardabellas (1963), who summarized data from a large
number of lysimeters:
R = 1.1∙P − ΔP

(10)

where P is annual precipitation and ΔP depends on soil and land cover.
For sandy soils with vegetation, ΔP = 433 mm, and the corresponding
results are shown in Table 4. Eq. (10) provides recharge estimatessignificantly larger than the numerical simulations and WTF applied toour site, except for the driest year, for which the agreement is reason­
able. It should be noted that R in Eq. (10) corresponds to the drainage
flux measured in lysimeters, which can be significantly larger thanrecharge, if there is a possibility of upward capillary flow or water up­
take by plant roots from depths below the bottom of the lysimeter.
Moreover, Eq. (10) does not account for the presence of low perme­ability layers, which was the case on our site.
Another simple formula relating recharge to precipitation overannual or longer time scales is:
R = ω∙P(11)

where ω = R/P is the recharge–precipitation ratio, which depends on
various factors. Pazdro & Kozerski (1990) suggested ω = 0.25 to 0.3 for
12

A. Gumuła-Kawęcka et al.

Journal of Hydrology 605 (2022) 127283

Fig. 8. Measured and simulated positions of the groundwater table (h, measured with respect to the profile bottom), representing the best fit for each profile obtainedusing simulations with lateral outflow – forest profiles P1 and P4.

glacial and fluvioglacial sands, such as those occurring on the Brda
outwash plain. The R/P ratio for the whole 3-year observation period
obtained using numerical simulations (0.24–0.28) and from WTF-MRC(0.20–0.30) is close to this estimate. In the more recent study by Duda
et al. (2011), ω was assumed to be a product of four partial coefficients:
ω = αβ γδ, representing the soil type (α), land cover (β), ground slope (γ),
and depth to groundwater table (δ). The suggested α value for sandy
deposits is 0.14 to 0.22, while β is 1 for grassland and 0.9 for forests. In
our case, the factors γ and δ were equal to 1, because the terrain isrelatively flat and the depth to the groundwater table is more than 2 m.
This gives a range of ω values from 0.14 to 0.22 for grassland and 0.13 to0.20 for the forest. The resulting recharge estimates are quite close to thenumerical simulations and WTF-MRC in the second (dry) and third(average) years, while they are significantly lower in the first (wet) year
and, consequently, also in terms of the total 3-year recharge. The co­efficients suggested by Duda et al. (2011) were inferred from the base­flow analysis for several Polish rivers. Regional (watershed) scale
recharge estimated from baseflow tends to be smaller than local-scaleestimations in areas with the groundwater table below the root zone

because a part of recharged groundwater is lost to evapotranspiration inriver valleys and wetlands before discharging to the river. However,recharge estimates for the region of our experimental site (an upper Brda
river watershed, Hobot et al., 2012) obtained from baseflow measure­ments are about 155 to 193 mm/yr (26% to 33% of the average pre­cipitation in Chojnice), which is in good agreement with the results ofour numerical simulations (an average of the 3-year period is 167 to 191
mm/yr).
Since the plant water demand does not change from one year to theother as much as precipitation, the fraction of total precipitation
reaching groundwater is smaller in dry years than in wet years(Kowalski, 2007). This is consistent with the results of numerical sim­
ulations, which indicate a much larger R/P ratio for the wettest year
(0.36–0.38) than the driest year (0.10–0.24). According to the WTFMRC method, the variability of the R/P ratio is smaller, and thelowest R/P ratio is in the third year, representing average conditions
between wet and dry years. Eq. (10) predicts a very strong dependenceof the R/P ratio on precipitation, which is not reflected in our results. In
contrast, Eq. (11) provides a constant ratio ω = R/P, i.e. the recharge is a13

A. Gumuła-Kawęcka et al.

Journal of Hydrology 605 (2022) 127283

Fig. 9. Measured and simulated positions of the groundwater table (h, measured with respect to the profile bottom), representing the best fit for each profile obtainedusing simulations with lateral outflow – grassland profiles P2 and P3.

fixed fraction of annual precipitation regardless of the actual amount ofprecipitation. Neither of these equations appeared adequate in our case.
However, Eq. (11) seems to work reasonably well in dry or average
years.
The increasing occurrence of extremely wet and dry periods alongwith temperature rise leading to higher evapotranspiration and shorter
snow season are the main symptoms of climate change in Central Europe(Goudie, 2006; IPCC, 2021), which in turn affect groundwater recharge(eg., Haidu and Nistor, 2019; Jaworska-Szulc, 2015; Neukum andAzzam, 2012;). Our modeling approach seems to be a promising tool toevaluate recharge under different historical or future climate change
scenarios (Gumuła-Kawęcka et al. 2021). In particular, it is possible torepresent groundwater table fluctuations using a 1D variably-saturated
flow model, if the parameters describing lateral outflow are known.
They can be obtained from limited time series of water table measure­ments, and the model can be further enhanced by including timevariable water level in the surface water body (lake or stream).

4.8. Influence of plant cover
It could be expected that recharge in a forest is smaller than in a
grassland area due to higher interception by the tree canopy and theability of tree roots to extract water from larger depths. However, there
does not appear to be a universal relationship between recharge inforests and grasslands, since it may depend on particular plant species.
According to Dripps & Bradbury (2010), recharge under the coniferous
forest was only slightly smaller than under grass cover. Grinevskii &Novoselova (2011) found that recharge under the forest was actuallylarger than under cropland.
In our simulations, recharge depended strongly on LAI. On average,
recharge estimates for the grassland profiles were larger than for theforest profiles (Fig. 6) because the assumed range of LAI values washigher for the forest than for grassland, despite a partial overlap.
However, looking only at the three scenarios, which showed the best
agreement with the observed water table fluctuations, we can see thatrecharge in the grassland profile P3 (LAI 2.9) was actually lower than in
P1 or P4 (with an LAI of 2.1). This is a consequence of our choice of LAI
14

A. Gumuła-Kawęcka et al.

Journal of Hydrology 605 (2022) 127283

Fig. 10. Sensitivity of the water table position in profile P1 (h, measured with respect to the profile bottom) to the parameters of lateral outflow boundary condition:
Hdr [m asl] (left) and C [(cm.d)− 1 (right).

Fig. 11. Master recession curves (MRC) representing the relationship between the groundwater table elevation (H [m a.s.l.]) and daily decrements of the water level(ΔHMRC [mm]).

15

A. Gumuła-Kawęcka et al.

Journal of Hydrology 605 (2022) 127283

Fig. 12. Increments of groundwater table caused by recharge in profile P3, according to the RISE and MRC methods.
Table 4
Recharge (R) and recharge/precipitation ratios (ω = R/P) obtained from numerical simulations WTF (MRC) and simple estimation methods.
Simulations
2017/18
2018/192019/20
total

MRC

Eq. (10)Eq. (11)

R [mm]

R/P [-]

R [mm]

R/P [-]

R [mm]

R/P [-]

R [mm]

R/P [-]

330–352
44–111
86–152
501–573

0.36–0.38
0.10–0.24
0.13–0.23
0.24–0.28

184–322
84–148
102–176
410–606

0.20–0.35
0.18–0.32
0.15–0.26
0.20–0.30

579
76
304
959

0.63
0.16
0.46
0.47

120–20260–102
87–147
267–452

0.13–0.22
0.13–0.22
0.13–0.22
0.13–0.22

values for the simulations, with the upper estimate for grass higher thanthe lower estimate for pine forest, as explained in Section 3.3. Our re­sults emphasize the need for accurate estimations of LAI. However,identification of actual LAI values in the study area and their possible
evolution in time was outside the scope of the present work.
In contrast, WTF did not show significant differences in recharge
between grassland and forest profiles. The results of the WTF method are
affected by the land cover in a larger area, not only in the close vicinityof the monitoring well. Forests are intermixed with open grasslands,
crop fields, and small built-up regions in our study area. Consequently,
water table fluctuations in both monitoring wells showed a similar
pattern, representing average recharge.

condition in numerical simulations, one can check if the vadose zone
model is consistent with the observed groundwater table fluctuations.
An additional benefit is the calculation of specific yield from the simu­
lation results. The specific yield can be later used in the WTF-MRC
analysis to see if the recharge estimates are consistent with those ob­tained from simulations.
One can note a conceptual similarity between the numerical simu­
lations with the modified bottom boundary condition and the masterrecession curve (MRC) approach. In both cases, the aim is to include
lateral groundwater flow in the analysis of point-scale observations. In
this way, we can obtain a more accurate description of the vadose zone –groundwater interactions without the need to solve more complex 2D or
3D models.
The incorporation of groundwater level observations in vadose zone
models significantly reduces the uncertainty in recharge estimations.
Our results showed that soil water measurements alone are insufficient
to constrain the models, especially if there is no detailed data on
evapotranspiration.
In view of the significant year-to-year variability of precipitation andrecharge, numerical simulations and WTF should preferably be applied
on longer time scales to achieve meaningful recharge estimations. In
numerical simulations, it is important to use a warm-up period with

5. Conclusion
This study highlights the need to use multiple methods and obser­
vation periods to increase confidence in recharge estimations, especiallygiven the increasing occurrence of extremely wet and dry periods intemperate climate zone. In our approach, numerical simulations ofvertical flow in the vadose zone were combined with the analysis ofwater table fluctuations. We showed how these two methods could
inform and cross-validate each other. Using a modified boundary
16

Journal of Hydrology 605 (2022) 127283

A. Gumuła-Kawęcka et al.

realistic weather data to accurately reproduce the initial state of soil
water.
Further research is warranted on the combined use of 1D numericalmodels and WTF. In particular, there is a need to investigate otheranalytical formulas describing the lateral outflow boundary condition
and to develop tools for calibrating their parameters jointly with thecalibration of soil hydraulic parameters and/or ET-related parameters.
We believe that the approach presented here can be used to improve
recharge estimates in unconfined aquifers over a wide range of
geological and hydrological settings.

De Silva, M.S., Nachabe, M.H., Šimůnek, J., Carnahan, R., 2008. Simulating root water
uptake from a heterogeneous vegetative cover. J. Irrig. Drain. Engin.-ASCE 134 (2),
167–174.
de Vries, J.J., Simmers, I., 2002. Groundwater recharge: an overview of processes andchallenges. Hydrogeol. J. 10 (1), 5–17.
Dripps, W.R., 2012. An integrated field assessment of groundwater recharge. OpenHydrol. J. 6 (1), 15–22.
Dripps, W.R., Bradbury, K.R., 2010. The spatial and temporal variability of groundwater
recharge in a forested basin in northern Wisconsin. Hydrol. Process. 24, 383–392.
Duda, R., Witczak, S., Żurek, A., 2011. Mapa wrażliwości wód podziemnych Polski na
zanieczyszczenie w skali 1:500 000 [Groundwater vulnerability map of Poland 1:500000]. Akademia Górniczo-Hutnicza im. S, Staszica w Krakowie, Kraków.
Dyck, P.P., Chardabellas, P., 1963. Wege zur Ermittlung der nutzbaren
Grundwasserreserven. Ber. Geol. Ges DDR 8, 245–262.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of Field Water Use and Crop
Yield. John Wiley & Sons, New York, NY.
Goudie, A.S., 2006. Global warming and fluvial geomorphology. Geomorphology 79 (34), 384–394.
Graf, R., Przybyłek, J., 2014. Estimation of shallow groundwater recharge using a GISbased distributed water balance model. Quaestiones Geographicae 33 (3), 7–37.
Grinevskii, S.O., Novoselova, M.V., 2011. Regularities in the formation of groundwater
infiltration recharge. Water Resour. 38 (2), 175–186.
Gumuła-Kawęcka, A., Jaworska-Szulc, B., Szymkiewicz, A., Gorczewska-Langner, W.,
2021. Impact of climate change on groundwater recharge in shallow sandy aquifer(Brda outwash plain, Pomeranian Region, Northern Poland). Subbmitted.
Gumuła-Kawęcka, A., Szymkiewicz, A., Jaworska-Szulc, B., Pruszkowska-Caceres, M.,
Gorczewska-Langner, W., Kaźmierczak, B., Kutyłowska, M., Piekarska, K.,
Jadwiszczak, P., 2018. Preliminary estimation of groundwater recharge on Brda
river outwash plain. E3S Web Conf. 44, 00050. https://doi.org/10.1051/e3sconf/
20184400050.
Haidu, I., Nistor, M.M., 2019. Long-term effect of climate change on groundwater
recharge in the Grand Est region of France. Meteorol. Appl. 27 (1), 1–17.
He, J., Hantush, M.M., Kalin, L., Rezaeianzadeh, M., Isik, S., 2021. A two-layer numericalmodel of soil moisture dynamics: model development. J. Hydrol. 602, 126797.
https://doi.org/10.1016/j.jhydrol.2021.126797.
Healy, R. W., 2010. Estimating groundwater recharge. Cambridge University Press,
Cambridge.
Heppner, C.S., Nimmo, J.R. 2005. A computer program for predicting recharge with amaster recession curve. U.S. Geological Survey Scientific Investigations Report 20055172.
Heppner, C.S., Nimmo, J.R., Folmar, G.J., Gburek, W.J., Risser, D.W., 2007. Multiplemethods investigation of recharge at a humid-region fractured Rock Site,
Pennsylvania, USA. Hydrogeol. J. 15 (5), 915–927.
Hobot, A., Banaszak, K., Stolarska, M., Sowińska, K., Serafin, R., Stachura, A., 2012.
Warunki korzystania z wód zlewni rzeki Brdy, Etap 1 – Dynamiczny bilans ilościowy
zasobów wodnych [Conditions of water usage in Brda river watershed, Stage 1 –Dynamic water budget], report for Regional Water Management Office in Gdańsk.
Pectore-Eco company, Gliwice.
Hölting, B., Coldewey, W.G., 2019. Hydrogeology. Springer.
Huet, M., Chesnaux, R., Boucher, M.-A., Poirier, C., 2016. Comparing various approaches
for assessing groundwater recharge at a regional scale in the Canadian Shield.
Hydrolog. Sci. J. 61 (12), 2267–2283.
Ines, A.V.M., Mohanty, B.P., 2008. Near-surface soil moisture assimilation for
quantifying effective soil hydraulic properties using genetic algorithm: 1. Conceptualmodeling. Water Resour. Res. 44, W06422.
IPCC, 2021. Summary for Policymakers. In: Climate Change 2021: The Physical Science
Basis. Contribution of Working Group I to the Sixth Assessment Report of theIntergovernmental Panel on Climate Change, Cambridge University Press. In Press.
Jackson, R.B., Canadell, J., Ehleringer, J.R., Mooney, H.A., Sala, O.E., Schulze, E.D.,
1996. A global analysis of root distributions for terrestrial biomes. Oecologia 108(3), 389–411.
Jagodziński, A.M., Kałucka, I., 2008. Age-related changes in leaf area index of young
Scots pine stands. Dendrobiology 59, 57–65.
Jaworska-Szulc, B., 2015. Impact of climate change on groundwater resources in a young
glacial multi-aquifer system. Pol. J. Environ. Stud. 24 (6), 2447–2457.
Jaworska-Szulc, B., Pruszkowska-Caceres, M., Szymkiewicz, A., Gumuła-Kawęcka, A.,
Gorczewska-Langner, W., 2017. Wahania zwierciadła płytkich wód podziemnych na
obszarze sandrowym w rejonie Chojnic [Groundwater table fluctuations in theshallow aquifer of the Brda River outwash plain]. Przegląd Geologiczny 65 (11),247–255.
Jie, Z., van Heyden, J., Bendel, D., Barthel, R., 2011. Combination of soil-water balance
models and water-table fluctuation methods for evaluation and improvement ofgroundwater recharge calculations. Hydrogeol. J. 19 (8), 1487–1502.
Johnson, A.I., 1967. Specific yield: compilation of specific yields for various materials.
US Geological Survey Water-Supply Paper 1662-D.
Kowalski, J., 2007. Hydrogeologia z podstawami geologii [Hydrogeology and elementsof geology], third ed. Wydawnictwo Uniwersytetu Przyrodniczego we Wrocławiu,
Wrocław.
Krogulec, E., 2010. Evaluation of infiltration rates within the Vistula Rivervalley, central
Poland. Acta Geol. Polonica 60 (4), 617–628.
Kurylyk, B.L., MacQuarrie, K.T.B., 2013. The uncertainty associated with estimating
future groundwater recharge: A summary of recent research and an example from asmall unconfined aquifer in a northern humid-continental climate. J. Hydrol. 492,
244–253.

CRediT authorship contribution statement
Anna Gumuła-Kawęcka: Conceptualization, Formal analysis,
Investigation, Methodology, Writing – original draft. Beata JaworskaSzulc: Conceptualization, Formal analysis, Investigation, Methodology.
Adam Szymkiewicz: Conceptualization, Formal analysis, Methodology,
Writing – original draft. Wioletta Gorczewska-Langner: Formal anal­ysis, Investigation. Małgorzata Pruszkowska-Caceres: Investigation.
Rafael Angulo-Jaramillo: Methodology, Investigation, Writing – re­view & editing. Jirka Šimůnek: Methodology, Software, Writing – re­view & editing.
Declaration of Competing Interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.
Acknowledgments
The work was funded by National Science Centre (NCN), Poland(grant 2015/17/B/ST10/03233 “Groundwater recharge on outwash
plain”). We highly appreciate the comments of four anonymous re­viewers, which helped to improve the original manuscript.
ReferencesÅberg, S.C., Korkka-Niemi, K., Rautio, A., Salonen, V.P., Åberg, A., 2019. Groundwater
recharge/discharge patterns and groundwater–surface water interactions in asedimentary aquifer along the River Kitinen in Sodankylä, northern Finland. Boreal
Environ. Res. 24, 155–187.
Ala-aho, P., Rossi, P.M., Kløve, B., 2014. Estimation of temporal and spatial variations ingroundwater recharge in unconfined sand aquifers using Scots pine inventories.
Hydrol. Earth Syst. Sci. Discuss. 11, 7773–7826.
Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998. Crop evapotranspiration-Guidelines
for computing crop water requirements-FAO Irrigation and drainage paper 56. Fao,
Rome, 300(9), D05109.
Batalha, M.S., Barbosa, M.C., Faybishenko, B., van Genuchten, M.T., 2018. Effect of
temporal averaging of 381 meteorological data on predictions of groundwater
recharge. J. Hydrol. Hydromech. 66, 143–152.
Batelaan, O., De Smedt, F., 2001. WetSpass: a flexible, GIS based, distributed recharge
methodology for regional groundwater modelling. IAHS PUBLICATION 11–18.
Beegum, S., Šimůnek, J., Szymkiewicz, A., Sudheer, K., Nambi, I., 2018. Updating thecoupling algorithm between HYDRUS and MODFLOW in the HYDRUS package for
MODFLOW. Vadose Zone J. 17 (1), 1–8.
Bouwer, H., 1989. Estimating and enhancing groundwater recharge. In: Sharma, M.L.
(Ed.), Grounwater recharge. A.A, Balkema, Rotterdam, pp. 1–10.
Braden, H., 1985. Ein Energiehaushalts- und Verdunstungsmodell for Wasser und
Stoffhaushaltsuntersuchungen landwirtschaftlich genutzer Einzugsgebiete.
Mittelungen Deutsche Bodenkundliche Geselschaft, 42, 294-299.
Brunner, P., Doherty, J., Simmons, C.T., 2012. Uncertainty assessment and implications
for data acquisition in support of integrated hydrologic models. Water Resour. Res.
48 (7), W07513.
Callahan, T.J., Vulava, V.M., Passarello, M.C., Garrett, C.G., 2012. Estimating
groundwater recharge in lowland watersheds. Hydrol. Process. 26 (19), 2845–2855.
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water
retention characteristics. Water Resour. Res. 24 (5), 755–769.
Cherkauer, D.S., Ansari, S.A., 2005. Estimating ground water recharge from topography,hydrogeology, and land cover. Ground Water 43 (1), 102–112.
Delin, G.N., Healy, R.W., Lorenz, D.L., Nimmo, J.R., 2007. Comparison of local-to
regional-scale estimates of ground-water recharge in Minnesota, USA. J. Hydrol. 334(1-2), 231–249.

17

A. Gumuła-Kawęcka et al.

Journal of Hydrology 605 (2022) 127283
and Multiple Solutes in Variably-Saturated Media. Version 4.17. Department ofEnvironmental Sciences. University of California Riverside, Riverside, CA, USA.
Šimůnek, J., van Genuchten, M., Šejna, M., 2008. Development and applications of the
HYDRUS and STANMOD software packages and related codes. Vadose Zone J. 7 (2),587–600.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2016. Recent developments andapplications of the HYDRUS computer software packages. Vadose Zone J. 15 (7).
Smerdon, B.D., Mendoza, C.A., Devito, K.J., 2008. Influence of subhumid climate andwater table depth on groundwater recharge in shallow outwash aquifers. Water
Resour. Res. 44, W08427.
Sophocleous, M.A., 1991. Combing the soil-water balance and water-level fluctuation
methods to estimate natural groundwater recharge: practical aspects. J. Hydrol. 124,
229–241.
Sousa, M.R., Jones, J.P., Frind, E.O., Rudolph, D.L., 2013. A simple method to assess
unsaturated zone time lag in the travel time from ground surface to receptor.
J. Contam. Hydrol. 144, 138–151.
Szilagyi, J., Zlotnik, V., Gates, J.B., Jozsa, J., 2011. Mapping mean annual groundwater
recharge in the Nebraska Sand Hills. USA. Hydrogeol. J. 19, 1503–1513.
Szymkiewicz, A., Gumuła-Kawęcka, A., Šimůnek, J., Leterme, B., Beegum, S., JaworskaSzulc, B., Pruszkowska-Caceres, M., Gorczewska-Langner, W., Angulo-Jaramillo, R.,
Jacques, D., 2018a. Simulation of the freshwater lens recharge using the HYDRUS
and SWI2 packages for MODFLOW. J. Hydrol. Hydromech. 66, 246–256.
Szymkiewicz, A., Gumuła-Kawęcka, A., Potrykus, D., Jaworska-Szulc, B., PruszkowskaCaceres, M., Gorczewska-Langner, W., 2018b. Estimation of conservative
contaminant travel time through vadose zone based on transient and steady flow
approaches. Water 10 (10), 1417.
Szymkiewicz, A., Savard, J., Jaworska-Szulc, B., 2019. Numerical analysis of recharge
rates and contaminant travel time in layered unsaturated soils. Water 11 (3), 545.
Szymkiewicz, A., Potrykus, D., Jaworska-Szulc, B., Gumuła-Kawęcka, A., PruszkowskaCaceres, M., Dzierzbicka-Głowacka, L., 2020. Evaluation of the influence of farming
practices and land use on groundwater resources in a coastal multi-aquifer system in
Puck region (northern Poland). Water 12 (4), 1042.
Twarakavi, N., Šimůnek, J., Seo, S., 2008. Evaluating Interactions between groundwaterand vadose zone using the HYDRUSbased flow package for MODFLOW. Vadose Zone
J. 7, 757–768.
van Dam, J.C., Huygen, J., Wesseling, J.G., Feddes, R.A., Kabat, P., van Walsum, P.E.V.,
Groenendijk, P., van Diepen, C.A., 1997. Theory of SWAP version 2.0, Report 71,
Dept. of Water Resour. Wageningen Agricultural University, Wageningen, theNetherlands.
van Genuchten, M.T., 1980. A closed form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.
Von Hoyningen-Hüne, J., 1983. Die Interception des Niederschlags inlandwirtschaftlichen Beständen. Schriftenreihe des DVWK 57, 1-53.
Wendland, F., Kunkel, R., Müller, U., Tetzlaff, B., 2008. Area – differentiated modellingof groundwater recharge rate for determining quantitative status of groundwater;
Case Study: Federal State of Lower Saxony, Germany. In: Mikulič, Z., Andjelov, M.
(Eds.) Groundwater modelling. Proceedings of invited lecturers of Symposium on
groundwater flow and transport modelling, Ljubljana, Slovenia, 28-31 January
2008, 73-79.
Wu, L., 2014. Groundwater recharge modelling approach to identify climate change
impacts using groundwater levels from Tärnsjö. TRITA-LWR Degree Project 14, 19.
Zendehboudi, S., Shafiei, A., Chatzis, I., Dusseault, M.B., 2012. Numerical simulation offree fall and controlled gravity drainage processes in porous media. J. Porous Media
15 (3), 211–232.
Zomlot, Z., Verbeiren, B., Huysmans, M., Batelaan, O., 2015. Spatial distribution ofgroundwater recharge and base flow: assessment of controlling factors. J. Hydrol.
Reg. Stud. 4, 349–368.

Lerner, D.N., Isaar, A.S., Simmers, I., 1990. Groundwater Recharge: A Guide toUnderstanding and Estimating Natural Recharge, International Contributions toHydrogeology. Verlag Heinz Heise, Hanover.
Leterme, B., Mallants, D., 2011. Climate and land use change impacts on groundwater
recharge. Proceedings ModelCARE2011. Models - Repositories of Knowledge. IAHS
Publ., Leipzig, Germany, 18-22 September 2011.
Leterme, B., Mallants, D., Jacques, D., 2012. Sensitivity of groundwater recharge using
climatic analogues and HYDRUS-1D. Hydrol. Earth Syst. Sci. 16 (8), 2485–2497.
Liu, Y., Yamanaka, T., Zhou, X., Tian, F., Ma, W., 2014. Combined use of tracer approachand numerical simulation to estimate groundwater recharge in an alluvial aquifer
system: a case study of Nasunogahara area, central Japan. J. Hydrol. 519, 833–847.
Muter, K., 2002. Mapa Hydrogeologiczna Polski w skali 1:50 000, ark. Cekcyn (204)[Hydrogeological Map of Poland 1:50 000, sheet Cekcyn (204)]. PIG, Warszawa.
Nastev, M., Rivera, A., Lefebvre, R., Martel, R., Savard, M., 2005. Numerical simulationof groundwater flow in regional rock aquifers, southwestern Quebec, Canada.
Hydrogeol. J. 13 (5-6), 835–848.
Neukum, C., Azzam, R., 2012. Impact of climate change on groundwater recharge in asmall catchment in the Black Forest, Germany. Hydrogeol. J. 20 (3), 547–560.
Nimmo, J.R., Horowitz, Ch., Mitchell, L., 2015. Discrete-storm water-table fluctuation
method to estimate episodic recharge. Groundwater 53(2), 282–292.
Niswonger, R., Prudic, D., Regan, R., 2006. Documentation of the Unsaturated-Zone Flow(UZF1) Package for modeling unsaturated flow between the land surface and thewater table with MODFLOW-2005: U.S. Geological Techniques and Methods Book 6,
Chapter A19.
Oficjalska, H., Gregosiewicz, R., 2000. Mapa Hydrogeologiczna Polski w skali 1:50 000,
ark. Tuchola (203) [Hydrogeological Map of Poland 1:50 000, sheet Tuchola (203)].
PIG, Warszawa.
Pazdro, Z., Kozerski, B., 1990. Hydrogeologia ogólna [Hydrogeology]. Wydawnictwa
Geologiczne, Warszawa.
Pozdniakov, S., Vasilevskiy, P., Grinevskiy, S., 2015. Estimation of groundwater recharge
by flow in vadose zone simulation at the watershed with different landscapes andsoil profiles. Engin. Geol. Hydrogeol. Bulgarian Acad. Sci. 29, 47–58.
Puchalski, T., Prusinkiewicz, Z., 1990. Ekologiczne podstawy siedliskoznawstwa leśnego[Ecological foundations of forest habitats]. PWRiL, Warszawa.
Rutkowski, P., Samborska, K., Wajsowicz, T., Konatowska, M., Budka, A., MaciejewskaRutkowska, I., Rybarczyk, A., 2017. Rozmieszczenie korzeni trzech gatunków drzew
w glebach Słowińskiego Parku Narodowego [Distribution of roots of three tree
species in soils of Słowiński National Park]. Acta Scientiarum Polonorum. Silvarum
Colendarum Ratio et Industria Lignaria 16 (4).
Sainju, U.M., Good, R.E., 1993. Vertical root distribution in relation to soil properties in
New Jersey Pinelands forests. Plant Soil 150 (1), 87–97.
Scanlon, B.R., Healy, R.W., Cook, P.G., 2002a. Choosing appropriate techniques forquatifying groundwater recharge. Hydrogeol. J. 10, 18–39.
Scanlon, B.R., Christman, M., Reedy, R.C., Porro, I., Simunek, J., J., Flerchinger, G.N.,
2002b. Intercode comparisons for simulating water balance of surficial sediments insemiarid regions. Water Resour. Res. 38 (12), 1323.
Schroeder, P.R., Dozier, T.S., Zappi, P.A., McEnroe, B.M., Sjostrom, J.W., Peyton, R.L.,
1994. Hydrologic Evaluation of Landfill Performance (HELP) model. Engineering
Documentation for Version 3.0. U.S. Environmental Protection Agency, Risk
Reduction Engineering Laboratory, Cincinnati, OH, USA.
Scurlock, J. M. O., Asner, G. P., & Gower, S. T. (2001). Worldwide historical estimates ofleaf area index, 1932–2000. ORNL/TM-2001/268, 34.
Şen, Z., 2015. Applied Drought Modelling, Prediction, and Mitigation. The, first ed.
Elsevier, Amsterdam.
Simmers, I., 1997. Recharge of Phreatic Aquifers in (Semi-) Arid Areas. A. A Balkema,
Rotterdam.
Šimůnek, J., Šejna, M., Saito, H., Sakai, M., van Genuchten, M.T., 2013. The HYDRUS-1D
Software Package for Simulating theOne-Dimensional Movement of Water, Heat,18