Science of the Total Environment 877 (2023) 162904

Contents lists available at ScienceDirect

Science of the Total Environment
journal homepage: www.elsevier.com/locate/scitotenv

Impact of climate change on groundwater recharge in shallow young glacial
aquifers in northern Poland⁎

Anna Gumuła-Kawęcka a, , Beata Jaworska-Szulc a, Adam Szymkiewicz a, Wioletta Gorczewska-Langner a,
Rafael Angulo-Jaramillo b, Jirka Šimůnek cabcGdańsk University of Technology, Faculty of Civil and Environmental Engineering, ul. Gabriela Narutowicza 11/12, 80-233 Gdańsk, Poland
Univ Lyon, Université Claude Bernard Lyon 1, CNRS, ENTPE, UMR 5023 LEHNA, F-69518 Vaulx-enVelin, France
University of California Riverside, CA, USA

H I G H L I G H T S

G R A P H I C A L

A B S T R A C T• Signiﬁcant air temperature rise occurredin last 70 years in northern Poland.
• Climate change causes disturbed precipitation pattern and intensive rain events.
• Frequent rainy/dry years result in observed groundwater table decrease in thelast years.
• Decreasing recharge trends in 2000–2020
were conﬁrmed by numerical simulations.
• Tracer experiments induce signiﬁcantwater content inﬂuence on vadosezone ﬂow.

A R T I C L E

I N F O

Editor: Christian Herrera
Keywords:
Unsaturated zone
Groundwater recharge
Climate action
Numerical model
HYDRUS-1DA B S T R A C T

We investigated the inﬂuence of climate change in the period 1951–2020 on shallow aquifers in the Brda and Wda outwash plains (Pomeranian Region, Northern Poland). There was a signiﬁcant temperature rise (0.3 °C/10 years), which accelerated after 1980 (0.66 °C/10 years). Precipitation became increasingly irregular – extremely rainy years occurred rightafter or before extremely dry years, and intensive rainfall events became more frequent after 2000. The groundwater level
decreased over the last 20 years, even though the average annual precipitation was higher than in the previous 50 years.
We carried out numerical simulations of water ﬂow in representative soil proﬁles for the years 1970–2020 using the
HYDRUS-1D model, developed and calibrated during our earlier work at an experimental site in the Brda outwash
plain (Gumuła-Kawęcka et al., 2022). We used a relationship between the water head and ﬂux at the bottom of the soil
proﬁles (the third-type boundary condition) to reproduce groundwater table ﬂuctuations caused by recharge variability
in time. The calculated daily recharge showed a decreasing linear trend for the last 20 years (0.05–0.06 mm d−1/
10 years), and dropping trends in water table level and soil water content in the entire proﬁle of vadose zone. Field tracer
experiments were performed to estimate impact of extremely rain events on water ﬂux in vadose zone. The results suggestthat tracer travel times are strongly determined by water content in the unsaturated zone which is determined by precipitation amount in span of weeks, rather than extremely high precipitation events.

1. Introduction⁎ Corresponding author.
E-mail address: annkawec@pg.edu.pl (A. Gumuła-Kawęcka).

Accurate groundwater recharge estimation is essential in sustainable
water management to protect the quantity and quality of water resources.
This knowledge is required to quantify the safe yield of an aquifer in the

http://dx.doi.org/10.1016/j.scitotenv.2023.162904
Received 11 November 2022; Received in revised form 24 February 2023; Accepted 12 March 2023
Available online 17 March 2023
0048-9697/© 2023 Elsevier B.V. All rights reserved.

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

Science of the Total Environment 877 (2023) 162904

of change in the precipitation totals is above 10 mm y−1/10 years for the
winter and 5–10 mm y−1/10 years for the spring, which are the highest
values in Poland (Mager et al., 2009). On the contrary, rainfalls in the summer decreased, but the number of irregular, extreme rainfall events during
the dry periods increased (Czarnecka and Nidzgorska-Lencewicz, 2012;
IPCC, 2021; Mager et al., 2009). The continental inﬂuence is shrinking
and retreating because the climate is getting more maritime characteristics (Degimendžić et al., 2004; Żmudzka, 2009). In the last years, thenumber of cloudy days in the central Pomerania region has become
the largest in Poland (Żmudzka, 2007), which is caused by a change in
the atmospheric circulation, which determines the inﬂuence of marine
or continental air masses (Twardosz et al., 2011). As a result, a large variation in the annual precipitation totals is observed, and extremely rainy
years more often occur right before or right after extremely dry years(Somorowska, 2022b).
Recently, Gumuła-Kawęcka et al. (2018, 2022), Gumuła-Kawęcka
(2023) conducted studies on groundwater recharge in Cekcyn, near
Tuchola, located in the southern part of the Brda sandur (Fig. 1). Numericalsimulations of water ﬂow in the unsaturated zone using literature-based
soil hydraulic parameters showed signiﬁcant variability in annual recharge
in the period 2003–2016, which corresponded to year-to-year differencesin precipitation (Gumuła-Kawęcka et al., 2018). In the follow-up work,
Gumuła-Kawęcka et al. (2022) used a combination of water content measurements in four soil proﬁles, numerical simulations of vadose zone
ﬂow, and the water table ﬂuctuation analysis to estimate groundwater recharge in the period 2017–2020, which included an extremely wet year
followed by an extremely dry year. Estimates of annual groundwater recharge differed by a factor from 1.8 to 7.5 between the wet and dry years,
depending on the soil proﬁle and the recharge estimation method. Furthermore, joint use of water content and groundwater table position measurements allowed us to ﬁnd soil hydraulic parameters, which were the mostrepresentative at a larger spatial scale. Gumuła-Kawęcka et al. (2022) developed an approach that allowed reproducing groundwater table ﬂuctuations using a 1D model of vertical ﬂow in a soil proﬁle. This was achievedby extending the solution domain into the aquifer and using a relationshipbetween water head and ﬂux (the third-type boundary condition) to represent lateral groundwater ﬂow from the soil proﬁle towards nearby discharge zones (e.g., lakes or streams).
The present work builds on the results of our previous studies and aims
to investigate long-term trends in climate variables and their inﬂuence on
shallow aquifers in the Brda and Wda outwash plains. Our objectives are(i) to assess changes in temperature, precipitation, and potential evapotranspiration from 1951 to 2020 based on available climate data, (ii) to analyze historical (1951–1970) and recent (2002−2020) measurements ofgroundwater table ﬂuctuations, (iii) to perform water ﬂow simulations in
representative soil proﬁles for the period 1970–2020, to validate themodeling approach developed in our earlier work for long-term simulations, (iv) to evaluate the trends in groundwater recharge, water table ﬂuctuations, and near-surface soil water contents obtained from simulations,and (v) to investigate extreme rain events inﬂuence on vadose zone ﬂow
using tracer experiments.

hydrological balance of catchment (Bouwer, 1989; de Vries and Simmers,
2002; Sophocleous, 1991; Xiang et al., 2019). The groundwater recharge
ﬂux also plays a signiﬁcant role in the transport of contaminants because
it determines the travel time through the unsaturated zone (e.g., Beroet al., 2016; Selker et al., 1996; Sousa et al., 2013; Szymkiewicz et al.,2018, 2020), directions of groundwater ﬂow in shallow aquifers(e.g., Scanlon et al., 2002a) and particular groundwater intakes tohealth-related hazards. In the last decades, climate change has affected
groundwater resources in many parts of the world (e.g., Gong et al.,2012; Guevara-Ochoa et al., 2020; Hu et al., 2019; Hughes et al., 2021;
Meixner et al., 2016; Qiu et al., 2023; Rodríguez-Huerta et al., 2019;
Tauﬁk et al., 2019; Yifru et al., 2021), and Central Europe is not an exception (e.g., Eckhardt and Ulbrich, 2003; Hattermann et al., 2008; JaworskaSzulc, 2015; Neukum and Azzam, 2012; van Roosmalen et al., 2007;
Wegehenkel and Kersebaum, 2009; Woldeamlak et al., 2007). Shallow
sandy aquifers, including aquifers located in outwash plains, are particularly vulnerable because they are recharged mainly by rainfall and snowmelt inﬁltration. An outwash plain, also named ‘sandur’, is a sandy
deposit formed by glacial ablation water. It is a geomorphological form
widely occurring in the glacial landscape of the North European Plain,
Fennoscandia, northern Syberia, Canada, and the northern US. A signiﬁcantspatial extent and high hydraulic conductivity of the outwash deposits
make sandurs abundant groundwater bodies commonly used to supply
water. The presented investigation was conducted in the Brda river sandur(Pomeranian region, Northern Poland), which is connected with the adjacent outwash plain of the Wda river and forms the largest system of shallow
aquifers in Poland.
Climate change in Central Europe is mainly related to the temperaturerise and an increase in evapotranspiration, changes in the amount and seasonal distribution of precipitation, and a shorter snow season (Goudie,
2006). Since 1980's there has been a noticeable increase in average temperature in Poland (Degimendžić et al., 2004; Kożuchowski et al., 2000;
Kundzewicz and Radziejewski, 2002; Mager et al., 2009; Żmudzka,2009). The latest IPCC report shows a signiﬁcant temperature rise ofabout 0.5 °C per decade in northern Poland from 1980 to 2015 (Gutiérrez
et al., 2021). The trend is not uniform for the entire territory of Poland. According to Mager et al. (2009), the Brda outwash plain area is where themost intensive climate warming occurs – above 0.25 °C/10 years in the period 1966–2006. The temperature increase is the highest in winter(0.35 °C–0.45 °C per 10 years) and spring (0.25–0.35 °C/10 years). The temperature rise in summer and autumn is more moderate (0.15–0.25 °C/
10 years). Rising temperatures cause an increase in potential evapotranspiration (PET) (Bogawski and Bednorz, 2016; Nistor et al., 2022;
Somorowska, 2022a; Teuling et al., 2009), which in turn contributes to adecrease in soil moisture content, a reduction of groundwater recharge,
and lower groundwater levels. Several signiﬁcant droughts occurred inEurope in the last years (e.g., Bastos et al., 2020; Christian et al., 2020;
García-Herrera et al., 2010; García-Herrera et al., 2019; Grillakis, 2019;
Ionita et al., 2017). Especially, high temperatures during the growing season and snow cover reductions resulted in severe soil moisture deﬁcits inCentral Europe (Boergens et al., 2020; Hari et al., 2020; Ionita et al.,2020; Ladányi et al., 2021; Pińskwar et al., 2020; Somorowska, 2022b).
The pattern of precipitation change is more complicated. The linear
trend for northern Poland in the period 1980–2015 published in the IPCC
report varies in a wide range from 0 mm y−1/10 years to an increase ofabout 44 mm y−1/10 years (Gutiérrez et al., 2021). It is difﬁcult to identify
a general trend in the rainfall change in Poland because it is different ineach region (Czarnecka and Nidzgorska-Lencewicz, 2012; Kożuchowski,
2004; Mager et al., 2009; Ziernicka-Wojtaszek, 2006; Żmudzka, 2002).
One of the highest precipitation rises is observed in the northern Poland –in some regions exceeding 10 mm y−1/10 years in the period 1966–2006(Mager et al., 2009). In the last years, the precipitation amount during
the winter increased in entire northern Poland (Czarnecka andNidzgorska-Lencewicz, 2012; Mager et al., 2009). The rise of minimal temperature results in replacing snowfalls with rainfalls, which changes theprecipitation pattern (Wibig and Glowicki, 2002). The average coefﬁcient

2. Materials and methods
2.1. Study area - Brda and Wda outwash plains
The Brda and Wda river sandurs were formed during the Weichselian(Vistulian/Wisconsinan) glaciation (115,000–11,700 years ago) as wide
outwash plains by meltwater emerging from the glacial lobe in the north(Fig. 1). They were drained by a system of braided streams and rivers,
which gradually transform into a single-channel river in the ToruńEberswalde ice-marginal valley in the south (Mojski, 2005). The deposits
consist of sand supported by gravel and pebbles. Fragments of glacial till,
eluted from the adjacent moraines, occur above the sandur deposits or between the outwash layers. The outwash sediments are mostly underlain
by glacial till or older outwash material (sand and gravel).
2

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

Science of the Total Environment 877 (2023) 162904

Fig. 1. Study area – the Brda and Wda outwash plains (modiﬁed from Gumuła-Kawęcka et al., 2018).

Daily data on maximum and minimum temperature, air humidity, windspeed, and cloud cover are available for the period 1966–2020. This data
allows for calculating the PET from the Penman-Monteith formula, assuming the Leaf Area Index is equal to 2.10, considered a representative value
for the area in our earlier study (Gumuła-Kawęcka et al., 2022).
Monthly and yearly average temperatures were evaluated using a statistical method introduced by Miętus et al. (2002) in climate change investigations in the northern part of the Pomeranian region, including the Wda
outwash plain. The method is based on an empirical percentile rank,which was applied to classify long-term data from the Chojnice station(years 1951–2020) to select outstanding values. Eleven percentile ranks
were distinguished to use as a thermal classiﬁcation (Table 1, Suppl.
Tables 1–2). Monthly and yearly rainfall totals were analyzed using themethod proposed by Miętus et al. (2005), similar to the one used for temperatures (Table 1, Suppl. Tables 3-4).
Daily precipitation were also analyzed to determine occurrence of extremely rain events for decades. Due to data availability records from1971 and 2020 was used. Percentile rank for entire 50-years period was calculated to ﬁnd days with precipitation exceeding 90th (4.90 mm/d), 95th
(8.20 mm/d) and 99th percentile (18.64 mm/d). We estimated part ofintensive rainfalls in overall precipitation total and number of days withextremely abundant rain for each decade.

The climate of the entire Poland territory is determined by maritime air
masses from the west, a polar climate of the Arctic, continental air from theeast, and subtropical fronts from southern Europe and the Atlantic Ocean.
The average air temperature in the Brda sandur region is 7.4 °C and the annual amount of precipitation is 581 mm/year.1 The climate demonstrates
an inﬂuence of the Baltic Sea with chilly springs (6.8 °C), cool summers(16.5 °C), and warm autumns (7.9 °C). The winters are moderately cold (averaging −1,6 °C) with poor snow cover. The average precipitation total
during winter months is 104 mm/year. The most intensive rainfalls occur
in the summer when the average sum of rainwater reaches 217 mm. The
precipitation total is higher in autumn (on average 140 mm) than in spring(on average 119 mm).
2.2. Climate data
Climate data were obtained from the Chojnice meteorological station(N 53.717, E 17.550) (Fig. 1), operated by the Polish Institute of Meteorology and Water Management (IMGW-PIB). Only monthly data (average temperature and precipitation totals) are available for the period 1951–1965.
1
Based on measurements from the Chojnice meteorological station (IMGW-PIB) during theyears 1951–2020.

3

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

Science of the Total Environment 877 (2023) 162904

Table 1
Annual average temperatures [°C] and precipitation totals [mm/y] at the Chojnice weather station classiﬁed based on percentile ranks.

measurements from 8 monitoring wells situated in the Brda sandur
and the adjacent Wda sandur (Fig. 1), published in the HydrologicalAtlas of Poland (Stachy, 1987). We chose wells not affected by groundwater extraction or ﬂuctuations in surface water bodies so that thewater level changes could be associated primarily with recharge variability.

2.3. Groundwater table data
Given the lack of water table measurements for the entire 70-year
period, the groundwater levels in the Brda and Wda sandurs were analyzed in two time periods: 1951–1970 and 2002–2019. The analysis
of groundwater levels in 1951–1970 was conducted based on
4

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

Science of the Total Environment 877 (2023) 162904

using daily meteorological data from Chojnice. We used two soil proﬁlesfrom our earlier study (Gumuła-Kawęcka et al., 2022), representing the forest (P1) and grassland (P2) areas of the experimental ﬁeld near Cekcyn on
the Brda outwash plain (Fig. 1). Detailed characterization of these proﬁles,including ﬁeld measurements, and set-up and calibration of the numericalmodels, is described in Gumuła-Kawęcka et al. (2022), and only a summaryis provided here. The shallow groundwater in entire experimental site isrecharged mainly by precipitation which was conﬁrmed with ﬁeld investigations (water level measurements and chemical analysis of groundwaterand surface waters) presented in Gumuła-Kawęcka (2023). In this study results of long-term simulations in the proﬁles were used to approximate recharge estimation for the whole experimental site divided into two zones
representing (I) forest – 1,164,062 m2, and (II) grassland – 206,250 m2(Fig. 1).
The HYDRUS-1D computer program (Šimůnek et al., 2013) was used tosolve the Richards equation describing vertical ﬂow in each of the soil proﬁles:Table 2
Correlation between precipitation, simulated groundwater recharge in P1 – P2, potential evapotranspiration (PET) and temperature based on monthly data. Correlation
scale bar: dark red – very strong/perfect (>0.9), light red – strong (0.7–0.9), orange –moderate (0.4–0.7), yellow – weak (0.2–0.4), green – no association (<0.2).
Precipitaon

Recharge
1P

2P

PET, LAI=2.1

Temperature-0,314-0,141

Precipitaon

1,000

0,048

0,439

1P

0,048

1,000

0,832

0,189

0,195

2P

0,439

0,832

1,000

-0,002

0,079

PET, LAI=2.1-0,3140,189-0,002

1,000

0,814

Temperature-0,141

0,195

0,079

0,814

1,000

Recharge

Table 3Time lag between the precipitation and the recharge in proﬁles P1 (woodland) andP2 (grassland) obtained with cross correlation between monthly precipitation andrecharge totals.
Time
periodProﬁleMean annual
effective
precipitation[mm/y]

Mean
annualrecharge[mm/y]

Time lag between
the precipitationand the recharge[months]

Correlation
coefﬁcient

1971–1975

P1
P2
P1
P2
P1
P2
P1
P2
P1
P2
P1
P2
P1
P2
P1
P2
P1
P2
P1
P2

108
141
207
226
122
146
127
156
143
176
154
177
205
240
223
253
98
136
178
217

97
170
94
193
98
161
81
159
67
168
93
176
130
224
115
224
107
177
104
211

17
3
9
2
6
5
24
5
19
5
7
5
6
3
7
3
6
4
7
5

0.26
0.32
0.38
0.50
0.17
0.30
0.31
0.40
0.40
0.49
0.47
0.55
0.38
0.48
0.35
0.47
0.23
0.43
0.37
0.52

1976–1980
1981–1985
1986–1990
1991–1995
1996–2000
2001–2005
2006–20102011–2015
2016–2020



∂θðhÞ
∂
∂h
∂k ðhÞ
¼k ðhÞþ
 SðhÞ
∂t
∂z
∂z
∂z

(1)

where θ is the volumetric water content [L3L−3], t is time [T], h is the waterpressure head (negative in the unsaturated zone), z is the spatial coordinate[L], k(h) is the hydraulic conductivity function of the unsaturated medium[LT−1], and S(h) is a sink function representing water uptake by plant roots[L3L−3 T−1].
Each soil proﬁle is 800 cm deep, extends below the groundwater table,
and consists of several soil layers. Hydraulic characteristics of each soil material were described with the van Genuchten-Mualem model (vanGenuchten, 1980):
θ−θr−m¼ ½1 þ ðαjhjÞng  g
θs −θrh  i2
1=m mgk ¼ ks kr ¼ ks Se0:5 1− 1−Se g
Se ¼

ð2aÞ
ð2bÞ

where Se is the effective saturation [−], θr is the residual water content[L3L−3], θs is the saturated water content [L3L−3], α is the parameter related to the average pore size [L−1], ng and mg are parameters related tothe pore size distribution [−], mg = 1–1/ng, ks is the vertical saturated hydraulic conductivity [L T−1], and kr is the relative hydraulic conductivity[−]. The saturated and residual water contents were estimated from ﬁeld
measurements and literature data according to Gumuła-Kawęcka et al.
(2022) (Suppl. Table 5). Other parameters of the van Genuchten functions(α, ng, and Ks) presented in Suppl. Table 5 were determined using inverse
modeling by ﬁtting the numerical solution to measured water contents ineach proﬁle at four depths (20, 50, 150, and 450 cm).
At the bottom of each proﬁle, we used the third-type boundary condition to represent lateral groundwater ﬂow towards a nearby surface water
body, provided by the “ﬂow to drains” option in HYDRUS-1D. We applied
the following relationship between the groundwater depth in the proﬁleand the horizontal ﬂux:

Fluctuations of groundwater levels in 2002–2020 were investigated inobservation wells of the Polish Geological Institute – National Research Institute (PGI–PIB) monitoring network located in shallow unconﬁned sandy
aquifers in the Brda and Wda outwash plains (Fig. 1). Monthly measurements from wells 183, 222, and 351/5 were published in the
Hydrogeological Annual Reports since November 2002. A preliminary
analysis of groundwater levels in the Brda sandur for the years
2002–2016 was published in Jaworska-Szulc et al. (2017).
2.4. Numerical simulations
Numerical simulations of groundwater recharge were performed for the51 years 1970–2020, with warm-up simulations for the years 1966–1969,

qdr ¼

4K h
ðH  H dr Þ2
L2dr

(3)

Table 4
Results of tracer experiments in proﬁles P1 (forest) and P2 (grassland).
Date ofinjection

ProﬁlenoSensor depth(z) [m]Time of penetration todepth z (t)Average water content invadose zone proﬁle [−]Precipitation totalin time t [mm]Number of extremely rain
events (intensity [mm/d])

July 3, 2017
Nov. 3, 2017
July 3, 2017
Nov. 3, 2017
April 27, 2018
Nov. 8, 2018

1P
1P
2P
2P
2P
2P

6.17
6.53
1.50
1.50
1.50
1.50

382
645
28
88
320
332

0.18
0.15
0.22
0.22
0.16
0.18

890.6
942.2
288.6
283.6
434.4
619.2

10 (18.4–46.6)
7 (18.4–46.6)
4 (19.8–30.0)
1 (18.4)
3 (19.2–46.6)
7 (19.2–46.8)

5

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

Science of the Total Environment 877 (2023) 162904

where qdr is the water ﬂux caused by ﬂow to drains [LT−1], Kh is the horizontal hydraulic conductivity of the aquifer [LT−1], Ldr is the horizontal distance from the soil proﬁle to the drain [L], H is the water table elevation in
the soil proﬁle [L], and Hdr is the water table elevation in drains [L], assumed equal to the elevation of the drain bottom. While the above equation
is typically used to describe ﬂow to drains, it represents a general formulafor horizontal groundwater ﬂow based on Dupuit's assumption.
In our earlier work (Gumuła-Kawęcka et al., 2022), we showed thatusing Eq. (3) as the bottom boundary condition allowed us to reproduce
natural ﬂuctuations of the groundwater table using the 1D vertical ﬂow
model in the soil proﬁle spanning the vadose zone and the upper part of
the aquifer. The parameters Hdr and C = 4Kh/(Ldr)2 are, in principle, unknown and must be calibrated using measurements of the groundwatertable level. In this study, we used the parameters calibrated earlier for our
experimental site: Hdr = 730 cm below the soil surface and C =1.32∙10−5 cm−1 d−1. Using a constant value of Hdr in time is a simpliﬁcation, as the surface water level in the receiving surface water body (groundwater discharge zone) is likely to vary during the simulation period.
However, such a simpliﬁcation can be justiﬁed, as no data on surface
water ﬂuctuations are available, and our analysis focuses on recharge andgroundwater table position trends rather than their speciﬁc values.
An atmospheric boundary condition was speciﬁed on the soil surface
based on daily weather data from the Chojnice station (daily maximum
and minimum temperatures, sunshine hours, air humidity, wind speed,
and precipitation totals). No water ponding was allowed (instantaneous
runoff), and the minimum pressure head on the soil surface in dry periods
was set to −1000 m. Potential evapotranspiration was calculated using
the Penman-Monteith equation (Allen, 1986). The leaf area index (LAI)
was assigned equal to 2.1 in the pine forest and 1.71 in the grassland, according to Gumuła-Kawęcka et al. (2022), using values from the range reported by Scurlock et al. (2001) for similar vegetation. The interception
was included according to the model implemented in HYDRUS-1D, basedon Braden (1985), van Dam et al. (1997), and Von Hoyningen-Hüne(1983) (Šimůnek et al., 2013). Root water uptake of grass and pine forest
was estimated using the Feddes et al. (1978) macroscopic model. The
plant-speciﬁc stress response function parameters were taken from De
Silva et al. (2008). In the absence of site-speciﬁc investigation of the root
zone, we assumed nonlinear root distribution with depth, following
Jackson et al. (1996), and the depth of the root zone was assigned according to Meyer et al. (1997) as 1.5 m for pine forest and 0.5 m for grass cover.
Since the initial conditions were not known, we started the simulations
with an arbitrary initial condition (a hydrostatic pressure head distribution
from the bottom of the proﬁle up to 2 m above the average position of thewater table and a uniform pressure head h = −2 m in the upper part of the
proﬁle) and used a warm-up period of four years 1966–1969 to minimize
the inﬂuence of the initial condition on the actual simulation period(1970–2020).
Snow accumulation and melting were included in the simulations. Thisrequired solving the heat transport equation in addition to the water ﬂow
equation. We used daily air temperatures as the surface boundary condition
for the heat equation, while at the bottom of the proﬁle, the temperature
gradient was set to zero. Soil thermal conductivity coefﬁcients were speciﬁed according to the Chung and Horton (1987) model, with parametersof the respective soil layers set to the default values for sand and loam, as
implemented in Hydrus 1D. A detailed description of the heat transportmodel and snow hydrology can be found in the Hydrus 1D manual(Šimůnek et al., 2013).

using a water can. The vegetative cover on both sites (pine forest in P1
and grass in P2) was left untouched in order to include its inﬂuence inthe experiments. The movement of chloride ions was tracked by measuring
the electrical conductivity (EC) using TDR sensors placed at depth 20, 60,
150 and 450 cm below the ground level (P2 only) and baro-divers in piezometers (P1 and P2) (Fig. 1). In both proﬁles the water content was recorded hourly with FDR sensors (P1) and TDR sensors (P2) at depth 20,
60, 150 and 450 cm for entire time of experiment. Sensors were placed inwell-permeable sandy layers (depth 20, 60 and 450 cm) and in semipermeable sandy loam (150 cm). In the same time precipitation totals
were collected by weather station located near proﬁle P2.
3. Results and discussion
3.1. Observed climate change trends
Statistically signiﬁcant increasing linear trends were ﬁtted to both temperature and PET data. Annual average temperatures in 1951–2020 corresponds to an increase of 0.3 °C/10 years (Suppl. Fig. 1), which is consistent
with Mager et al. (2009). The temperature rise for the last years is more intensive: 0.66 °C/10 years for the period 1980–2020, and 0.70 °C/10 years after
2000. This value exceeds the trend published by IPCC (Gutiérrez et al., 2021),
which is caused by abnormally warm years 2018–2020, which are not included in the IPCC report. The temperature increase results in proportionallyhigher potential evapotranspiration (Fig. 2, Suppl. Fig. 1). There is a linear
rise of PET during the 55 years (1966–2020) of about 25 mm y−1/
10 years. However, an increase in PET is about 33 mm y−1/10 years from1980 to 2020 and about 35 mm y−1/10 years for the period after 2000.
Two time periods with different annual temperatures can be distinguished(Table 1): the period between 1951 and 1987 (37 years), which was predominantly cold, and the period 1988–2020 (33 years), which was apparently
warmer. In the ﬁrst period, 25 years were classiﬁed as cold, including
6 years abnormally or extremely cold and 9 years very cold. Only 4 years
were warm, and two of them were slightly warm. On the contrary, in the second period (1988–2020), 22 years were evaluated as warm, and merely
3 years were cold. Six years were classiﬁed as abnormally or extremely
warm years, and 5 of them occurred in the last 7 years (2014–2015,
2018–2020). In both distinguished periods (1951–1987 and 1988–2020),
the number of years classiﬁed as “average” is relatively low (8 in each period).
Averages and standard deviations calculated for 5-year periods (Suppl.
Table 6) agree with the above ﬁndings. The temperature rise has been continuous for the last 30 years. The average 5-year temperatures before 1986
vary between 6.3 and 7.4 °C with a standard deviation (σ) between 0.4 and
0.9 °C. Since the late 1980s, the temperature rise has been constant andreached the maximum in the last ﬁve years (9.1 °C). The standard deviation
had similar values except for 1986–1990, when 1989 was abnormallywarm, and 1996–2000, when 2000 was extremely warm.
Analysis of annual precipitation totals in 1951–2020 recorded at the
Chojnice station show no statistically signiﬁcant trend. It is consistent
with the IPCC reports (Gutiérrez et al., 2021) which present wide rangeof potential annual rainfalls change. A numer of studies conducted for
Poland pointed that changes in the precipitation pattern are complex anddiffer depending on the region (Czarnecka and Nidzgorska-Lencewicz,
2012; Kożuchowski, 2004; Mager et al., 2009; Ziernicka-Wojtaszek, 2006;Żmudzka, 2002). Increasing intensive storms in the summer and replacingsnow with rain in winter bring difﬁculties to the long-term precipitation evaluation. Statistical analysis of daily rainfall totals between 1970 and 2020shows that precipitation are more intensive after 2000 (Fig. 2). The rain
events exceeding 90th percentile in 2001–2020 take 14.3 %–16.4 % of thetotal precipitation which is an apparent rise comparing to period
1971–2000 (11.3 %–12.8 %). The part of the outstanding rainfalls, above
99th percentile is also higher after 2000 (3.5 %–5.0 %) than in previous
years (2.4 %–3.5 %). Extremely rain events were also more frequent. In
2001–2020 387–436 d/10y bring precipitation exceeding 90th percentile, including 36–55 d/10y with rainfalls above 99th percentile. For 1971–2000
they were 318–380 d/10y and 26–36 d/10y, respectively. The results are

2.5. Tracer experiments
Chloride ions were applied as a tracer in the vadose zone because of
their conservative character and low background concentrations. The
tracer was applied in four injections around proﬁle P1 and P2 (July 3,
2017; November 3, 2017; April 27, 2018 and November 8, 2018) (Fig. 1).
In each case, 10 dm3/m2 of the chloride solution (25 g/dm3 NaCl) was
spilled on the ground surface over the area of 5 m2 (P1) or 7 m2 (P2)6

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

Science of the Total Environment 877 (2023) 162904

Fig. 2. Extremely precipitation totals (% in overall precipitation) and numer of days with heavy rainfalls (% per number of days with precipitation >0.1 mm/d) for each
decade between 1971 and 2020.

70 years, and abnormally wet/dry years were rare. As a result, the watertable position should be less disturbed by weather conditions. By comparison, in 2002–2019, climatic ﬂuctuations were signiﬁcant, including extraordinary low or high rainfalls observed along with the temperature rise.

consistent with Zeder and Fischer (2020) which reported more intensive
short- and long-duration extreme precipitation events in Germany,
Netherlands, Switzerland and part of Austria due to climate change.
Beside the daily data analysis, monthly and annual precipitation totals
were investigated using method by Miętus et al. (2005). Five classes weredistinguished based on the percentile rank of the long-term data from theyears 1951–2020 (Table 1, Suppl. Table 3). In a 70-year period, 26 years
were classiﬁed as average, 22 years were evaluated as dry, including 7 extremely dry years, and an equal number of rainy years (22) was observed,
including 7 extremely rainy years. By coincidence, the numbers of rainy
and dry years are equal (22), similar to the numbers of extremely rainy
and extremely dry years. Interestingly, the number of abnormal years (44wet or dry years) exceeds the number of average years (26).
The 5-year averages of yearly precipitation totals during the ﬁrst two
decades (1951–1970) vary between 527.0 and 623.3 mm/y (Suppl.
Table 6). The standard deviation was relatively low in 1951–1965(75.0–84.9 mm) and slightly higher in 1966–1970 (118.6 mm), caused
by abundant rainfalls in 1967 and 1970. After 1970, the occurrence of extremely wet and extremely dry years increased, and there are several periods when extremely wet years occur right after or right before dry or
extremely dry years (e.g., 1970–1972, 1979–1982, 2001–2003,
2006–2007, and 2014–2018). The standard deviation was signiﬁcant during 1976–1985 (141.9–146.7 mm) and 2001–2010 (about 135 mm). In
the last years, rainfall ﬂuctuations were particularly high; σ for
2014–2018 reached 156 mm (not shown in Suppl. Table 6). Although
long-term precipitation totals (Table 1, Suppl. Fig. 1) do not demonstrate
a clear trend, the abnormal years occurred more frequently in the years
1970–1990 (mostly dry years) and between 2000 and 2020 (mostly rainy
years). This is probably the reason why 5-year averages of yearly rainfalls
are slightly higher after 2000 (536.6–708.8 mm/y) than in the previous
years (489.8–623.3 mm/y) (Suppl. Table 6).

3.2.1. 1951–1970
The average annual depth of the water table in each year was comparedwith annual rainfall totals obtained from the Chojnice meteorological station (Suppl. Fig. 2). The results show that in all cases, the groundwaterlevel increased between 1951 and 1970. In general, the water table was
the lowest in 1954–1955, probably due to poor precipitation in 1951
(507 mm/y) and 1953 (425 mm/y). Even though in 1952, rainfall wasabove average (634 mm/y), it was not enough to replenish the aquifers to
the previous level. In the following years, precipitation was more abundant,with the annual rainfall exceeding 600 mm/y in 1956, 1958, 1960–1963,
1966–1967, and 1970. That caused a signiﬁcant increase in groundwaterlevels in all the wells. The highest water tables were observed in
1961–1963 and 1967–1968. The water table ﬂuctuations in the Brda
outwash plain in 1951–1970 were extraordinarily high: 463 cm inMąkowarsko, 168 cm in Osielsko, 152 cm in Pruszcz Pomorski, and104 cm in Sierosław. In the wells located in the Wda sandur (Czarna
Woda, Śliwice), an increase in groundwater levels was much lower(26–55 cm) during that time. Between 1951 and 1970 annual precipitation
totals were about average – the mean was 579 mm/y. It caused a rising
trend in the water table level in all monitoring wells because rainfalls/
snowfalls were relatively similar and regular. The presented analysis
shows that, even though extremely wet/dry years were rare during the analyzed period (1951–1970), water table differences in some cases reached>1.5 m, meaning the Brda outwash plain is highly vulnerable to climatic
ﬂuctuations.
3.2.2. 2002–2020
The analysis of groundwater level ﬂuctuations in the last two decades
shows that the water table response to precipitation takes from less thanhalf a year in shallow aquifers to about 1.5 years in deeper aquifers(Jaworska-Szulc et al., 2017). Apparent water table increases were

3.2. Observed groundwater table ﬂuctuations
In 1951–1970, the annual average temperature was relatively stable,
differences between annual precipitation totals were the lowest in the7

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

Science of the Total Environment 877 (2023) 162904

value. In P2, recharge variability is less pronounced, with a max/min
ratio about 4. This is still larger than the ratio between the maximum(835 mm/y) and minimum (312 mm/y) annual precipitation in the simulation period, which is 2.7. Range of simulated actual evapotranspiration iseven lower, between 309 and 533 mm/y in P1 and 285–512 mm/y in P1,
which gave maximum/minimum ratio of about 1.7–1.8. Higher variabilityin P1 can be attributed to larger interception and root depth in the forest, as
well as the lower permeability of the sandy loam layer, compared to P2.
These factors signiﬁcantly decreased water inﬁltration and percolation, especially in drier periods (e.g., Acharya et al., 2017; Awada et al., 2013;
Eldridge et al., 2015; Wine et al., 2015).
The ratio of average recharge to average precipitation (wR) is 0.17 in P1
and 0.32 in P2, which is consistent with values obtained in the literature forshallow sandy aquifers in a temperate climate (e.g., Dripps and Bradbury,
2010; Graf and Przybyłek, 2014; Krogulec, 2010; Leterme et al., 2012;
Pozdniakov et al., 2015; Wendland et al., 2008). However, the rechargeprecipitation ratios calculated for each year show strong variability, ranging from 0.04 to 0.47 in P1 and from 0.17 to 0.61 in P2 (again, variability
is larger in the forest proﬁle). For further analysis of dependence between
recharge, temperature, PET and precipitation a statistical correlation of annual data was carried (Table 2). Strong dependence was found between recharge simulated in proﬁles P1–P2, and between temperature and PET. The
relations between recharge(P1)–precipitation, recharge(P1)–temperatureand recharge(P1)–PET are weak or none. The correlation between recharge(P2)–precipitation is moderate, while for recharge(P2)–temperatureand recharge(P2)–PET there is no association. Presented analysis showsthat dependence between recharge and climatic factors (precipitation, temperature and PET) is not straightforward and recharge probably depends
also on meteorological conditions from previous years. This dependency
is expected to increase with increasing water table depth and decreasing
soil permeability (eg., Cook et al., 2003; Rossman et al., 2014; Wu et al.,1996).
The pattern of temporal changes in annual recharge generally mimics
the changes in precipitation, while the maxima and minima are shifted intime, especially in the P1 proﬁle (Fig. 3). The results show that the lowest
annual recharge does not necessarily occur in years classiﬁed as extremely
dry (Suppl. Tables 7–8). In proﬁle P1, 6 out of 10 lowest values of annual
recharge occurred in normal or even rainy years. However, in most cases,
these years were preceded by dry years. In proﬁle P2, the minimum recharge was more aligned with dry years, but 4 out of 10 lowest values
still corresponded to years classiﬁed as normal, rainy or even extremely
rainy. The larger discrepancy observed in P1 is related to a longer water
movement through the proﬁle. Similar observations can be made about
the maximum annual recharge values. A signiﬁcant proportion of high annual recharge values corresponds to normal or even dry years, preceded byobserved after the rainy summers in 2004, 2007, 2010, 2016, and 2017. Between autumn 2002 and spring 2007, the water table dropped, with only asmall increase in winter/spring 2005 caused by higher precipitation insummer 2004 (Suppl. Fig. 3). Abundant rainfalls in the summer of 2007 resulted in a temporary rise of shallow groundwater in 2008, after which thewater table gradually decreased in winter/spring 2010. Heavy rainfalls insummer/autumn 2010 caused signiﬁcant inﬂow to shallow aquifers between autumn 2010 and spring 2011. After that, the groundwater level
gradually decreased until it reached its lowest point in 2015/2016. Precipitation in summer 2012 slightly recharged the aquifers in winter/spring2013. Extraordinary high rainfalls in 2016–2017 (1403 mm in 2 years)
caused a rapid, even though temporary, water table rise in 2017–2018.
In the last two decades, a long-term decrease in groundwater levels has
been observed in the Brda outwash plain. This trend is apparent in monitoring wells observed for at least ten years in the entire Brda and Wda sandurs,
even though the annual precipitation totals are generally higher after 2000
(averagely 632 mm/y) than in the previous 50 years. Given the large distances between the analyzed piezometers and no signiﬁcant human impact
on groundwater intake in the nearby vicinity of each well, water table ﬂuctuations were probably caused by climate change (less snow, higher evapotranspiration, and more intensive rainfalls). Decreasing water table levels
were also observed in Germany (Hattermann et al., 2008; Natkhin et al.,2012; Zielhofer et al., 2022), Sweden and Finland (Nygren et al., 2020,2021). Similarly, a spring discharge reduction was observed in the
Odra basin (south-eastern Poland) (Tarka and Staśko, 2010), Bavaria(Germany), Lower Silesia (Poland), and northern Carpathians (Slovakia)(Staśko and Buczyński, 2018). The water table decrease is prognosed formost regions in Germany until the end of XXI century (Wunsch et al., 2022).
3.3. Numerical modeling results
3.3.1. Groundwater recharge
Annual groundwater recharge obtained in numerical simulations is presented in Fig. 3 and Suppl. Tables 7–8. The recharge was calculated as thewater ﬂux at the actual position of the water table, which varied in time.
The average annual recharge in the 51-year simulation period varies between 24 and 259 mm/y in proﬁle P1 (forest) and 94–364 mm/y in proﬁleP2 (grassland), with averages 99 mm/y and 189 mm/y, respectively(Fig. 3). For whole experimental site the recharge reaches between
27,588–301,492 m3/y in zone I (forest) and 19,470–74,972 m3/y in zone
II (grassland) (Suppl. Fig. 4). Obtained values are an approximate volumes
calculated based on recharge simulated in representative proﬁles P1–P2
and forest and grassland area estimated from topographic map. There is signiﬁcant year-to-year variability in annual recharge, especially in proﬁle P1,
where the maximum value is almost 11 times larger than the minimum

Fig. 3. Annual groundwater recharge in proﬁles P1–P2 obtained in numerical simulation compared to precipitation recorded between 1970 and 2020 (P = 586 mm/y on
average). Mean (RA), minimum (RMIN) and maximum recharge (RMAX) is expressed in mm/y and recharge/precipitation ratio (RA/P) has no units.
8A. Gumuła-Kawęcka et al.

Science of the Total Environment 877 (2023) 162904

Fig. 4. Simulated groundwater recharge trends in proﬁles P1 – P2 between 2000 and 2020. Decreasing recharge trends for the last two decades (2000–2020) are similar inboth cases – 0.05 mm d−1/10 years at P1 and 0.06 mm d−1/10 years at P2.

reached 8.9 °C, which demonstrated in higher mean evapotranspiration(412 mm/y for P2 and 446 mm/y for P1). However, the yearly precipitation (599 mm/y) was slightly above the long term average, which led tosimilar mean effective precipitation (152–187 mm/y).
The numerical simulation results show that despite the effective precipitation was similar for both periods, the average recharge was lower in
2014–2019 (89–179 mm/y) than in 1979–1984 (109–186 mm/y)(Fig. 5). The largest difference was observed in proﬁle P1, which was probably related to higher interception and evapotranspiration in the forest (P1)than in the grassland (P2). The recharge/precipitation ratio (wR) was from
0.15 to 0.30 lower in 2014–2019 than in 1979–1984. Therefore, eventhough the rainfall total was higher in 2014–2019, increased evapotranspiration reduced part of rainwater available for inﬁltration, which is demonstrated by lower wR.
The recharge time distribution is presented in Fig. 5. In general, verticalﬂow in proﬁle P2 is larger and more intensive than in proﬁle P1. It is probably caused by differences in the vegetation cover (mainly interception androot zone thickness) and more permeable sandy loam in proﬁle P2 (Suppl.
Table 5). In P1, groundwater recharge is lower, varies more smoothly overtime, and is slightly retarded compared to P2. In some cases, the verticalﬂow through the unsaturated zone can be very rapid, e.g., 17 mm/d in proﬁle P2 on July 17–18, 1980, after heavy rainfalls at the end of June 1980.
However, single intensive rain events, such as the storm on June 15,
2016, do not cause any outstanding effects on groundwater recharge. Although there are differences in vertical ﬂow through the unsaturated
zone in proﬁles P1 and P2, the time variability pattern of groundwater recharge is similar in both cases (Fig. 3–5).

rainy or extremely rainy years. Again, this effect is more pronounced in P1(7 out of 10 years) than in P2 (4 out of 10 years).
The delay in recharge arrival to the groundwater table with respect tothe time of water inﬁltration into soil depends on the depth of the watertable and the hydraulic conductivity and capacity of unsaturated soil.
Table 3 presents results of cross correlation between monthly effective precipitation (calculated as a difference between rainfall total and actual
evapotranspiration) and recharge in proﬁles P1 – P2 for years
1971–2020, analyzed in 5-years periods. In general, the time differences
are higher in proﬁle P1 (6–24 months) than in P2 (2–5 months). The largest
time differences are observed in the woodland proﬁle (P1), while in thegrassland (P2) they are more uniform. The analysis show that the time
lag do not depend directly on precipitation or recharge amount (Table 3)
but probably other factors such as, rainfall distribution and intensity,
water content in vadose zone or activation of preferential ﬂow channels
in the root zone are also important (e.g., Gee et al., 1994; Qiao et al.,2015, 2017; Tang, 1996; Wilcox, 2002).
Despite the high year-to-year recharge variability, some long-term tendencies can be found. Fig. 4 shows simulated daily recharge totals in proﬁles P1 – P2 between 2000 and 2020. Statistically signiﬁcant trends forthe last two decades is dropping in both proﬁles, with a decrease of
0.05 mm d−1/10 years at P1 and 0.06 mm d−1/10 years at P2. The falling
tendency in groundwater recharge agrees well with regional studies presented in Section 3.2.
A more detailed investigation of groundwater recharge was carried out
for two selected 6-year periods: 1979–1984 and 2014–2019 (Fig. 5). In
each of them, the annual precipitation was highly variable, from extremely
low to extremely high: 311.8–813.7 mm/y between 1979 and 1984, and433.4–818.9 mm/y between 2014 and 2019 (Table 1). During the ﬁrst period, starting in 1979 the mean annual rainfall (525 mm/y) was about
50 mm below the long-term average. The temperature was relatively low(7.0 °C on average) which resulted in mean evapotranspiration ranging
from 356 mm/y (P2) to 380 mm/y (P1). Thus, the average effective precipitation for simulated proﬁles was between 145 and 171 mm/y. The second
period (2014–2019) was extremely warm, the average temperature

3.3.2. Groundwater level
Daily groundwater level ﬂuctuations obtained for years 2000–2020
using numerical simulations are shown in Fig. 6. The variability patternin each proﬁle is consistent with recharge variations discussed earlier. Fluctuations in P1 are more smoothed and delayed in time compared to P2. As
mentioned above, this can be attributed to the slower water movement in
the P1 proﬁle. Statistically signiﬁcant linear trends ﬁtted to the data show9

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

Science of the Total Environment 877 (2023) 162904

Fig. 5. Daily groundwater recharge in proﬁles P1–P2 obtained using numerical simulations for 1979–1984 and 2014–2019. Precipitation (P), average (RA), minimum (RMIN)
and maximum recharge (RMAX) is expressed in mm/y and recharge/precipitation ratio (RA/P) has no units.

in vadose zone proﬁle at both the forest (P1) and grassland (P2) sites. The
decrease is larger in proﬁle P1 (0.0019–0.0098 per 10 years) than in P2(0.0006–0.0042 per 10 years). A decreasing trends in soil moisture is consistent with recent trends observed in the entire Poland (Pińskwar et al.,2020; Somorowska, 2022b).

decreasing tendency in both proﬁles, of about 43 mm/10 years in P1 and
33 mm/10 years in P2. This is consistent with water table measurementsin the recent period, described in Section 3.2. The results of numerical simulations should be considered qualitative rather than quantitative indicators due to the simpliﬁed representation of water table ﬂuctuations in the
1D vertical ﬂow model. Nevertheless, applying the third-type boundary
condition representing horizontal ﬂow in the aquifer seems to be a usefulapproach for preliminary assessment of the inﬂuence of climatic changes
on groundwater table positions in shallow aquifers, which are affected
mainly by natural hydrological processes.

3.4. Inﬂuence of extremely rain events on shallow groundwater
Two EC peaks at the water table level were observed in P1, while in the
proﬁle P2 apparent EC peaks were recorded by two sensors placed at thedepths 20 cm and 150 cm. Sensor at 50 cm below the ground level showed
irregular, increased values without a clear peak. Deeper sensors (450 cm
and at the water table) have not recorded any sign of the tracer. It is probably a result of limited ﬂow through the loamy sand layer at depth
200–320 cm, which is conﬁrmed by low values of water content at thedepth of 450 cm (averagely about 0.06).
The results of experiment show a signiﬁcant difference between years
2017 and 2018–2019 (Table 4). Depending on the moisture the tracer
travel time varied signiﬁcantly from 382 to 645 days in P1 and from28 days 332 days in P2. In 2017 the average water content in the vadosezone proﬁle was relatively high (0.18 in P1 and 0.22 in P2) and 4 extremelyrain events occurred during the experiment. The tracer during that time
moved through the vadose zone much faster than in the following period

3.3.3. Soil water content
Fig. 7 presents simulated soil water contents at depths of 20 cm, 50 cm,
150 cm and 450 cm for the period 2000–2020. Soil moisture contents in theforest (P1) ﬂuctuate between 0.025 and 0.19 in the top soil, 0.025–0.13 in
the upper sandy layer (Sand I), 0.25–0.41 in sandy loam and 0.04–0.07 in
the lower sandy layer (Sand II). In the grassland (P2) water contents oscillate between 0.025 and 0.24 in Sand I, 0.36–0.38 in sandy loam and
0.05–0.07 in Sand II. They vary signiﬁcantly throughout each year, being
highly inﬂuenced by dynamic conditions at the soil-atmosphere interface.
Statistically signiﬁcant linear trend lines were ﬁtted to each set of watercontent values sampled with a 1-day interval in the whole 21-year period.
The result shows that water contents are dropping at all measured depths
10

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

Science of the Total Environment 877 (2023) 162904

Fig. 6. Fluctuations of the groundwater table in proﬁles P1 and P2 obtained in numerical simulations.
Statistically signiﬁcant linear trends ﬁtted to the data show decreasing tendency in both proﬁles, of about 43 mm/10 years in P1 and 33 mm/10 years in P2.

(Table 4). The ﬁrst injection was applied 6 days after the intensive rainfall
in June 29 (38.8 mm/d) and the tracer in P2 moved to the depth 150 cm in
28 days. During the second experiment started on November 2017, only
one intensive precipitation was recorded. The tracer travel time was longer,
but still relatively short (88 days to depth 150 cm). In contrast, the third(April 2018) and fourth injection (November 2018) were applied during
a dry period, when water content was about 0.13 in P1 and 0.16–0.18 in
P2. In that case, it took 320–332 days for the tracer to reach the sensor at
depth 150 cm in P2, even though 3 extremely rain events were observed
during the third experiment, and 7 during the fourth experiment. Similartendency was noticed in proﬁle P1. The chlorides spilt in July 2017 appeared at the water table level after 382 days from the injection. Time
travel of the tracer introduced in November 2017 was much longer andreached 645 days to the water table at depth 653 cm.
Obtained tracer travel times under high water content are relatively
short, which can be also caused by some preferential ﬂow channels activated by heavy rainfalls. The largest differences in the results were obtained
for the shallow soil, up to 150 cm below the ground level (P2), which ismore affected by rainfalls changeability. Experiments conducted in southern Poland with lysimeters contained uncovered medium sand also gavehigh tracer velocities, from 35,7 mm/day to 15,4 mm/day (Żurek andMościcki, 2017). Similar incidents of rapid ﬂow in shallow unconﬁned
sandy aquifers were also described in other studies (Dripps, 2012;
Herrmann and Duncker, 2008; Seiler et al., 2002).
Presented results show that tracer travel times are strongly determined bywater content in the unsaturated zone which is controlled by precipitation
amount in span of weeks, rather than extremely high precipitation events.
Regular, several-days rainfalls increase the water content in entire vadosezone proﬁle, which improves percolation to the aquifer. On the contrary, single intensive rain events during dry periods inﬂuence mostly shallow subsurface part of unsaturated zone and has low importance to the groundwater
recharge. Decreasing trends of water content described in Section 3.3 may
cause a water ﬂux retardation and reduce the groundwater recharge.

3.5. Limitations and uncertainty
Our study shows impact of meteorological factors inﬂuenced by climate change, on unsaturated zone and shallow aquifer. We do not consider change in land use or groundwater resources uptake which wereminimal in our experimental ﬁeld. 1D numerical simulations performed
for representative vadose zone proﬁles, which results were extrapolated
for entire experimental ﬁeld to give an approximate recharge results inlateral scale. Some uncertainty is related to that approach, especially according to soil parameters estimations (α, ng, and Ks), LAI, root zone parameters and water level in adjacent lakes. However, recent ﬁeld
investigations in the experimental site with a variety of methods described in Gumuła-Kawęcka et al. (2022), Gumuła-Kawęcka (2023)
were used to minimalize the uncertainty. Calculations using Richards
equation, van Genuchten-Mualem model, horizontal ﬂux equation,
Penman-Monteith equation, Feddes model and Chung and Horton
model implemented in the numerical code also brings some simpliﬁcations. Nevertheless, HYDRUS-1D was successfully applied to groundwater recharge estimations (e.g., Batalha et al., 2018; Grinevskii and
Novoselova, 2011; Pozdniakov et al., 2015; Scanlon et al., 2002b;Šimůnek et al., 2016; Szymkiewicz et al., 2018; Twarakavi et al.,2008), including climate change impact (Leterme et al., 2012; Letermeand Mallants, 2011).
Using tracer methods in unsaturated zone also provides uncertainty.
Flow path of the tracer is usually assumed to be one-dimensional or piston
ﬂow (Walker, 1998), however variability of sediment layers and heterogeneity of soil properties can disturb the ﬂux transporting the tracer (Healy,2010). Hydrodynamic dispersion and diffusion also inﬂuence the tracer
movement according to advection-dispersion equation. Tracers in unsaturated zone represent a point-scale estimation of water ﬂow in vertical proﬁle which strongly depend on small-scale processes in the unsaturated
zone, such as preferential ﬂow in macropores or water ponding over lesspermeable layers.
11

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

Science of the Total Environment 877 (2023) 162904

Fig. 7. Simulated water contents at a depth of 20 cm, 50 cm, 150 cm and 450 cm in proﬁles P1 – P2.
The decrease is larger in proﬁle P1 (0.0019–0.0098 per 10 years) than in P2 (0.0006–0.0042 per 10 years).

12

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

Science of the Total Environment 877 (2023) 162904

4. Conclusions

Appendix A. Supplementary data

The analysis of historical weather records showed a considerable increase in temperatures in the study area, with an increase of 0.3 °C/
10 years in 1951–2020 and 0.66 °C/10 years between 1980 and 2020,
which exceeds the trend published by IPCC. The rising temperature results
in increasing evaporation of about 25 mm/10 years in 1966–2020 and
33 mm/10 years in 1980–2020. Annual precipitation amounts demonstrateincreasing irregularity, while daily records show more frequent extremelyrain events after 2000. Decrease in the groundwater table level observedin shallow aquifers in the Brda and Wda outwash plains suggests decreasing
recharge. Numerical simulations performed for two soil proﬁles with calibrated parameters conﬁrmed the trends of decreasing recharge, groundwater table elevation, and water contents in vadose zone over the last 20 years.
There is signiﬁcant year-to-year variability in annual groundwater recharge, with a ratio between maximum and minimum values of 3.9 in thegrassland proﬁle (P2) and 10.8 in the forest proﬁle (P1). However, the decreasing recharge trends for the last two decades (2000−2020) are similar
in both cases – 0.05 mm d−1/10 years at P1 and 0.06 mm d−1/10 years at
P2. Correlation between recharge and climatic factors (precipitation, temperature and PET) is not straightforward and recharge probably depends
also on meteorological conditions from previous years. Results of tracer experiments induce that time lag between effective precipitation and recharge do not depend directly on this two variables but also other factors
such as, rainfall distribution and intensity, water content in vadose zone
or activation of preferential ﬂow channels in the root zone. The soil moisture determine water ﬂux through unsaturated zone, thus the water content
decrease caused by climate change can cause a ﬂow retardation and reduce
the groundwater recharge. Our study shows that employing 1D model ofvertical ﬂow through the soil proﬁle provides a complete description of interplay between climatic factors (precipitation, evapotranspiration) and hydrologic conditions (water content in the unsaturated zone, groundwater
recharge, water table ﬂuctuations, horizontal ﬂow in the aquifer). The calculated water table ﬂuctuations were in quantitative agreement with observations, conﬁrming the utility of our modeling approach to predict thelong-term trends for shallow aquifers.

Supplementary data to this article can be found online at https://doi.
org/10.1016/j.scitotenv.2023.162904.
References
Acharya, B.S., Hao, Y., Ochsner, T.E., Zou, C.B., 2017. Woody plant encroachment alters soilhydrological properties and reduces downward ﬂux of water in tallgrass prairie. Plant
Soil 414, 379–391.
Allen, R., 1986. Penman for all seasons. Proc. ASCE, J. Irrigation and Drainage Eng. 112,
pp. 348–368.
Awada, T., El-Hage, R., Geha, M., Wedin, D.A., Huddle, J.A., Zhou, X., Msanne, J., Sudmeyer,
R.A., Martin, D.L., Brandle, J.R., 2013. Intra-annual variability and environmental controls over transpiration in a 58-year-old even-aged stand of invasive woody Juniperus
virginiana L. in the Nebraska Sandhills, USA. Ecohydrology 6, 731–740.
Bastos, A., Ciais, P., Friedlingstein, P., Sitch, S., Pongratz, J., Fan, L., Wigneron, J.P., Weber,U., Reichstein, M., Fu, Z., Anthoni, P., Arneth, A., Haverd, V., Jain, A.K., Joetzjer, E.,
Knauer, J., Lienert, S., Loughran, T., McGuire, P.S., Tian, H., Viovy, N., Zaehle, S.,
2020. Direct and seasonal legacy effects of the 2018 heat wave and drought on
European ecosystem productivity. Sci. Adv. 6 (2020), 1–14.
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.
Bero, N.J., Ruark, M.D., Lowery, B., 2016. Bromide and chloride tracer application to determine sufﬁciency of plot size and well depth placement to capture preferential ﬂow andsolute leaching. Geoderma 262, 94–100.
Boergens, E., Güntner, A., Dobslaw, H., Dahle, C., 2020. Quantifying the central European
droughts in 2018 and 2019 with GRACE follow-on. Geophys. Res. Lett. 47.
Bogawski, P., Bednorz, E., 2016. Atmospheric conditions controlling extreme summertime
evapotranspiration in Poland (central Europe). Nat. Hazards 81, 55–69.
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. MittelungenDeutsche Bodenkundliche Geselschaft. 42, pp. 294–299.
Christian, J.I., Basara, J.B., Hunt, E.D., Otkin, J.A., Xiao, X., 2020. Flash drought development
and cascading impacts associated with the 2010 Russian heatwave. Environ. Res. Lett. 15.
Chung, S.-O., Horton, R., 1987. Soil heat and water ﬂow with a partial surface mulch. Water
Resour. Res. 23 (12), 2175–2186.
Cook, P.G., Jolly, I.D., Walker, G.R., Robinson, N.I., 2003. From drainage to recharge to discharge: some timelags in subsurface hydrology. In: Alsharhan, A.S., Wood, W.W. (Eds.),
Developments in Water Science. vol. 50. Elsevier, pp. 319–326.
Czarnecka, M., Nidzgorska-Lencewicz, J., 2012. Wieloletnia zmienność opadów sezonowych
w Polsce. Woda-Środowisko-Obszary Wiejskie. 12, 2(38), pp. 45–60.
De Silva, M., Nachabe, M.H., Š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 and challenges. Hydrogeol. J. 10, 5–17.
Degimendžić, J., Kożuchowski, K., Żmudzka, E., 2004. Changes of air temperature and precipitation in Poland in the period 1951–2000 and their relationship to atmospheric circulation. Int. J. Climatol. 24, 291–310.
Dripps, W., 2012. An integrated ﬁeld assessment of groundwater recharge. Open Hydrol.J. 6,
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.
Eckhardt, K., Ulbrich, U., 2003. Potential impacts of climate change on groundwater recharge
and streamﬂow in a central European low mountain range. J. Hydrol. 284, 244–252.
Eldridge, D.J., Wang, L., Ruiz-Colmenero, M., 2015. Shrub encroachment alters the spatial
patterns of inﬁltration. Ecohydrology 8, 83–93.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of Field Water Use and Crop Yield.
John Wiley & Sons, New York, NY.
García-Herrera, R., Díaz, J., Trigo, R.M., Luterbacher, J., Fischer, E.M., 2010. A review of theEuropean summer heat wave of 2003. Crit. Rev. Environ. Sci. Technol. 40 (4), 267–306.
García-Herrera, R., Garrido-Perez, J.M., Barriopedro, D., Ordóñez, C., Vicente-Serrano, S.M.,
Nieto, R., Gimeno, L., Sorí, R., Yiou, P., 2019. The European 2016/17 drought. J. Clim.
32, 3169–3187.
Gee, G.W., Wierenga, P.J., Andraski, B.J., Young, M.H., Fayer, M.J., Rockhold, M.L., 1994.
Variations in water balance and recharge potential at three western desert sites. Soil
Sci. Soc. Am. J. 558 (1), 63–72.
Gong, J., Wang, K., Kellomäki, S., Zhang, C., Martikainen, P.J., Shurpali, N., 2012. Modelingwater table changes in boreal peatlands of Finland under changing climate conditions.
Ecol. Model. 244, 65–78.
Goudie, A., 2006. Global warming and ﬂuvial geomorphology. Geomorphology 79 (3–4),
384–394.
Graf, R., Przybyłek, J., 2014. Estimation of shallow groundwater recharge using a GIS-based
distributed water balance model. Quaestiones Geographicae 33 (3), 7–37.
Grillakis, M.G., 2019. Increase in severe and extreme soil moisture droughts for Europe underclimate change. Sci. Total Environ. 660, 1245–1255.
Grinevskii, S., Novoselova, M., 2011. Regularities in the formation of groundwater inﬁltration
recharge. Water Resour. 38 (2), 175–186.
Guevara-Ochoa, C., Medina-Sierra, A., Vives, L., 2020. Spatio-temporal effect of climatechange on water balance and interactions between groundwater and surface water inplains. Sci. Total Environ. 722, 137886.

CRediT authorship contribution statement
Anna Gumuła-Kawęcka: Conceptualization, Data curation, Formal
analysis, Investigation, Methodology, Software, Visualization, Writing –original draft. Beata Jaworska-Szulc: Conceptualization, Data curation,
Formal analysis, Investigation, Methodology, Supervision, Writing – review& editing. Adam Szymkiewicz: Conceptualization, Funding acquisition,
Formal analysis, Investigation, Methodology, Project administration, Supervision, Software, Writing – original draft, Writing – review & editing.
Wioletta Gorczewska-Langner: Data curation, Formal analysis, Methodology, Supervision, Visualization. Rafael Angulo-Jaramillo: Conceptualization, Writing – review & editing. Jirka Šimůnek: Software, Writing –review & editing.
Data availability
The authors do not have permission to share data.
Declaration of competing interest
The authors declare that they have no known competing ﬁnancial interests or personal relationships that could have appeared to inﬂuence thework reported in this paper.
Acknowledgments
This work has been supported by grant 2015/17/B/ST10/03233‘Groundwater recharge on outwash plain’ funded by National ScienceCentre (NCN), Poland.
13

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

Science of the Total Environment 877 (2023) 162904
Natkhin, M., Steidl, J., Dietrich, O., Dannowski, R., Lischeid, G., 2012. Differentiating between climate effects and forest growth dynamics effects on decreasing groundwater recharge in a lowland region in Northeast Germany. J. Hydrol. 448–449, 245–254.
Neukum, C., Azzam, R., 2012. Impact of climate change on groundwater recharge in a small
catchment in the Black Forest,Germany. Hydrogeol. J. 20, 547–560.
Nistor, M.-M., Satyanaga, A., Dezsi, .̧ S., Haidu, I., 2022. European grid dataset of actual evapotranspiration, water availability and effective precipitation. Atmosphere 13, 772.
Nygren, M., Giese, M., Kløveb, B., Haafa, E., Rossi, P.M., Barthel, R., 2020. Changes in seasonality of groundwater level ﬂuctuations in a temperate-cold climate transition zone.
J. Hydrol. X 8, 100062.
Nygren, M., Giese, M., Barthel, R., 2021. Recent trends in hydroclimate and groundwaterlevels in a region with seasonal frost cover. J. Hydrol. 602, 126732.
Pińskwar, I., Choryński, A., Kundzewicz, Z.W., 2020. Severe drought in the spring of 2020 in
Poland—more of the same? Agronomy 10, 1646.
Pozdniakov, S., Vasilevskiy, P., Grinevskiy, S., 2015. Estimation of groundwater recharge by
ﬂow in vadose zone simulation at the watershed with different landscapes and soil proﬁles. Eng.Geol.Hydrogeol. 29, 47–58.
Qiao, L., Zou, C.B., Will, R.E., Stebler, E., 2015. Calibration of swat model for woody plant encroachment using paired experimental watershed data. J. Hydrol. 523, 231–239.
Qiao, L., Zou, C.B., Stebler, E., Will, R.E., 2017. Woody plant encroachment reduces annual
runoff and shifts runoff mechanisms in the tallgrass prairie, USA. Water Resour. Res.
53, 4838–4849.
Qiu, J., Shen, Z., Xie, H., 2023. Drought impacts on hydrology and water quality under climate change. Sci. Total Environ. 858 (Part 1), 159854.
Rodríguez-Huerta, E., Rosas-Casals, M., Hernández-Terrones, L.M., 2019. Water societal metabolism in the Yucatan Peninsula. The impact of climate change on the recharge ofgroundwater by 2030. J. Clean. Prod. 235, 272–287.
Rossman, N.R., Zlotnik, V.A., Rowe, C.M., Szilagyi, J., 2014. Vadose zone lag time and potential 21st century climate change effects on spatially distributed groundwater recharge inthe semi-arid Nebraska Sand Hills. J. Hydrol. 519 (Part A), 656–669.
Scanlon, B.R., Healy, R.W., Cook, P.G., 2002a. Choosing appropriate techniques for quantifying groundwater recharge. Hydrogeol. J. 10, 18–39.
Scanlon, B.R., Healy, R.W., Cook, P.G., 2002b. Choosing appropriate techniques for quatifying
groundwater recharge. Hydrogeol. J. (10), 18–39.
Scurlock, J.M.O., Asner, G.P., Gower, S.T., 2001. Worldwide Historical Estimates of Leaf Area
Index, 1932–2000. ORNL/TM-2001/268, p. 34.
Seiler, K.-P., von Loewenstern, S., Schneider, S., 2002. Matrix and bypass-ﬂow in quaternary andtertiary sediments of agricultural areas in south Germany. Geoderma 105 (3–4), 299–306.
Selker, J.S., Steenhuis, T.S., Parlange, J.Y., 1996. An engineering approach to ﬁngered vadosepollutant transport. Geoderma 70 (2–4), 197–206.
Šimůnek, J., Šejna, M., Saito, H., Sakai, M., van Genuchten, M.Th, 2013. The HYDRUS-1D
Software Package for Simulating the One-Dimensional Movement of Water, Heat, andMultiple Solutes in Variably-Saturated Media. Version 4.17. Department of Environmental Sciences, University of California Riverside, Riverside, CA, USA.
Šimůnek, J., van Genuchten, M., Šejna, M., 2016. Recent developments and applications of
the HYDRUS computer software packages. Vadose Zone J. (15 (7)), 25.
Somorowska, U., 2022a. Changes in terrestrial evaporation across Poland over the past four
decades dominated by increases in summer months. Resources 11, 6.
Somorowska, U., 2022b. Ampliﬁed signals of soil moisture and evaporative stresses acrossPoland in the twenty-ﬁrst century. Sci. Total Environ. 812, 151465.
Sophocleous, M.A., 1991. Combing the soil-water balance and water-level ﬂuctuation
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.
Stachy, J., 1987. Atlas hydrologiczny Polski. Tom I. Instytut Meteorologii i Gospodarki
Wodnej, Wydawnictwa Geologiczne, Warszawa.
Staśko, S., Buczyński, S., 2018. Drought and its effects on spring discharge regimes in Poland
and Germany during the 2015 drought. Hydrol. Sci. J. (63(5)), 741–751.
Szymkiewicz, A., Gumuła-Kawęcka, A., Potrykus, D., Jaworska-Szulc, B., PruszkowskaCaceres, M., Gorczewska-Langner, W., 2018. Estimation of conservative contaminant
travel time through vadose zone based on transient and steady ﬂow approaches. Water
10 (10), 1417.
Szymkiewicz, A., Potrykus, D., Jaworska-Szulc, B., Gumuła-Kawęcka, A., PruszkowskaCaceres, M., Dzierzbicka-Głowacka, L., 2020. Evaluation of the inﬂuence of farming practices and land use on groundwater resources in a coastal multi-aquifer system in Puck region (northern Poland). Water 12 (4), 1042.
Tang, C., 1996. Interception and recharge processes beneath a Pinus Elliotii forest. Hydrol.
Process. 10 (11), 1427–1434.
Tarka, R., Staśko, S., 2010. Wahania zwierciadła wód podziemnych jako odzwierciedlenie
ekstremalnych sytuacji pogodowych. Wyjątkowe zdarzenia przyrodnicze na DolnymŚląsku i ich skutki…. Wydawnictwo Uniwersytetu Wrocławskiego, Wrocław,
pp. 241–255.
Tauﬁk, M., Veldhuizen, A., Wösten, J., van Lanen, H., 2019. Exploration of the importance ofphysical properties of Indonesian peatlands to assess critical groundwater table depths,
associated drought and ﬁre hazard. Geoderma 347, 160–169.
Teuling, A., Hirschi, M., Ohmura, A., Wild, M., Reichstein, M., et al., 2009. A regional perspective on trends in continental evaporation. Geophys. Res. Lett. 36 (2) American Geophysical Union.
Twarakavi, N., Šimůnek, J., Seo, S., 2008. Evaluating interactions between groundwater andvadose zone using the HYDRUSbased ﬂow package for MODFLOW. Vadose Zone J. 7,
757–768.
Twardosz, R., Niedźwiedź, T., Łupikasza, E., 2011. The inﬂuence of atmospheric circulationon the type of precipitation (Krakow, southern Poland). Theor.Apply Climatol. 104,
233–250.

Gumuła-Kawęcka, A., 2023. Groundwater Recharge Estimation on Outwash Plain Area Basedon Field Investigations and Numerical Modeling. Gdańsk University of Technology,Gdańsk, Poland Doctoral dissertation.
Gumuła-Kawęcka, A., Szymkiewicz, A., Jaworska-Szulc, B., Pruszkowska-Caceres, M.,
Gorczewska-Langner, W., 2018. Preliminary estimation of groundwater recharge on
Brda river outwash plain. 10th Conference on Interdisciplinary Problems in Environmental Protection and Engineering EKO-DOK 2018. EDP Sciences.
Gumuła-Kawęcka, A., Jaworska-Szulc, B., Szymkiewicz, A., Gorczewska-Langner, W.,
Pruszkowska-Caceres, M., Angulo-Jaramillo, R., Šimůnek, J., 2022. Estimation of groundwater recharge in a shallow sandy aquifer using unsaturated zone modeling and watertable ﬂuctuation method. J. Hydrol. 605, 127283.
Gutiérrez, J.M., Jones, R.G., Narisma, G.T., Alves, L.M., Amjad, M., Gorodetskaya, I.V., Grose,M., Klutse, N.A.B., Krakovska, S., Li, J., Martínez-Castro, D., Mearns, L.O., Mernild, S.H.,
Ngo-Duc, T., van den Hurk, B., Yoon, J.-H., 2021. Atlas. In: Masson-Delmotte, V., Zhai, P.,
Pirani, A., Connors, S.L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M.I.,
Huang, M., Leitzell, K., Lonnoy, E., Matthews, J.B.R., Maycock, T.K., Waterﬁeld, T.,
Yelekçi, O., Yu, R., Zhou, B. (Eds.), Climate Change 2021: The Physical Science Basis.
Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press.
Hari, V., Rakovec, O., Markonis, Y., Hanel, M., Kumar, R., 2020. Increased future occurrencesof the exceptional 2018–2019 Central European drought under global warming. Sci. Rep.
10, 12207.
Hattermann, F., Post, J., Krysanova, V., Conradt, T., Wechsung, F., 2008. Assessment of water
availability in a central-European river basin (Elbe) under climate change. Adv. Clim.
Change Res. 4, 1673–1719 (Suppl.).
Healy, R.W., 2010. Estimating Groundwater Recharge. Cambridge University Press.
Herrmann, A., Duncker, D., 2008. Runoff formation in a tile-drained agricultural basin of theHarz Mountain foreland,Northern Germany. Soil Water Res. 3, 83–97.
Hu, G., Zhao, L., Li, R., Wu, X., Wu, T., Xie, C., Zhu, X., Su, Y., 2019. Variations in soil temperature from 1980 to 2015 in permafrost regions on the Qinghai-Tibetan Plateau based on
observed and reanalysis products. Geoderma 337, 893–905.
Hughes, A., Mansour, M., Ward, R., Kieboom, N., Allen, S., Seccombe, D., Charlton, M.,
Prudhomme, C., 2021. The impact of climate change on groundwater recharge:
national-scale assessment for the British mainland. J. Hydrol. 598, 126336.
Ionita, M., Nagavciuc, V., Kumar, R., Rakovec, O., 2020. On the curious case of the recent decade, mid-spring precipitation deﬁcit in central Europe. npj Clim. Atmos. Sci. 3, 49.
Ionita, M., Tallaksen, L.M., Kingston, D.G., Stagge, J.H., Laaha, G., Van Lanen, Scholz, P.,
Chelcea, S.M., Haslinger, K., 2017. The European 2015 drought from a climatologicalperspective. Hydrol. Earth Syst. Sci. 21, 1397–1419.
IPCC, 2021. Summary for policymakers. Climate Change 2021: The Physical Science Basis.
Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental 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.
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 wód podziemnych na
obszarze sandrowym w rejonie Chojnic. Prz. Geol. 65 (11), 247–255.
Kożuchowski, K., 2004. Zmienność opadów atmosferycznych w Polsce w XX i XXI wieku. In:Kożuchowski, K. (Ed.), Skala, uwarunkowania i perspektywy współczesnych zmian
klimatycznych w Polsce. Biblioteka, Łódź, pp. 47–58.
Kożuchowski, K., Degimendžić, J., Fortuniak, K., Wybig, J., 2000. Trends to changes in seasonal aspects of the climate in Poland. Geographia Polonica. Polish Academy of Science,
pp. 7–24.
Krogulec, E., 2010. Evaluation of inﬁltration rates within the Vistula Rivervalley, central
Poland. Acta Geol. Pol. 60 (4), 617–628.
Kundzewicz, Z., Radziejewski, M., 2002. Detection of Changes of Climate and Hydrological
Processes (in Polish) Sorus, Poznań.
Ladányi, Z., Barta, K., Blanka, V., Pálffy, B., 2021. Assessing available water content of sandy
soils to support drought monitoring and agricultural water management. Water Resour.
Manag. 35, 869–880.
Leterme, B., Mallants, D., 2011. Climate and land use change impacts on groundwater recharge. Proceedings ModelCARE2011. Models - Repositories of Knowledge. IAHS Publ.,
Leipzig, Germany September 18-22.
Leterme, B., Mallants, D., Jacques, D., 2012. Sensitivity of groundwater recharge using climatic analogues and HYDRUS-1D. Hydrol. Earth Syst. Sci. 16, 2485–2497.
Mager, P., Kasprowicz, T., Farat, R., 2009. Changes in air temperature and precipitation in
Poland in 1966-2006. In: Leśny, J. (Ed.), Acta Agrophysica. Climate Change and Agriculture in Poland - Impact, Mitigation and Adaptation Measures. Polish Academy of Science,
Poznań, pp. 19–38.
Meixner, T., Manning, A.H., Stonestrom, D.A., Allen, D.M., Ajami, H., Blasch, K.W.,
Brookﬁeld, A.E., Castro, Ch.L., Clark, J.C., Gochis, D.J., Flint, A.L., Neff, K.L., Niraula,R., Rodell, M., Scanlon, B.R., Singha, K., Walvoord, M.A., 2016. Implications of projectedclimate change for groundwater recharge in the western United States. J. Hydrol. 534,
124–138.
Meyer, P., Rockhold, M., Gee, G., 1997. Uncertainty Analyses of Inﬁltration and Subsurface Flow and Transport for SDMP Sites. Paciﬁc Northwest National Laboratory,
Richland.
Miętus, M., Owczarek, M., Filipiak, J., 2002. Warunki termiczne na obszarze Wybrzeża i
Pomorza w świetle wybranych klasyﬁkacji. Materiały Badawcze IMGW, Seria
Meteorologia. 36, p. 56.
Miętus, M., Filipiak, J., Owczarek, M., Jakusik, E., 2005. Zmienność warunków opadowych
polskiego wybrzeża Morza Bałtyckiego w świetle kwantylowej klasyﬁkacji opadowej.
Materiały Badawcze IMGW, Seria Meteorologia 37, 59.
Mojski, J., 2005. Ziemie polskie w czwartorzędzie: Zarys morfogenezy. PIG, Warszawa.

14

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

Science of the Total Environment 877 (2023) 162904
Wunsch, A., Liesch, T., Broda, S., 2022. Deep learning shows declining groundwater levels inGermany until 2100 due to climate change. Nat. Commun. 13, 1221.
Xiang, W., Si, B.C., Biswas, A., Li, Z., 2019. Quantifying dual recharge mechanisms in deep unsaturated zone of Chinese Loess Plateau using stable isotopes. Geoderma 337, 773–781.
Yifru, B.A., Chung, I.-M., Kim, M.-G., Chang, S.W., 2021. Assessing the effect of land/use land
cover and climate change on water yield and groundwater recharge in East African Rift
Valley using integrated model. <sb:contribution><sb:title>J.</sb:title></sb:
contribution><sb:host><sb:issue><sb:series><sb:title>Hydrol. Reg. Stud.</sb:title></sb:series></sb:issue></sb:host> 37, 100926.
Zeder, J., Fischer, E.M., 2020. Observed extreme precipitation trends and scaling in Central
Europe. Weather Clim. Extremes 29, 100266.
Zielhofer, C., Schmidt, J., Reiche, N., Tautenhahn, M., Ballasus, H., Burkart, M., Linstädter, A.,
Dietze, E., Kaiser, K., Mehler, N., 2022. The lower Havel River region (Brandenburg,
Germany): a 230-year-long historical map record indicates a decrease in surface water
areas and groundwater levels. Water 14, 480.
Ziernicka-Wojtaszek, A., 2006. Zmienność opadów atmosferycznych na obszarze Polski w
latach 1971–2000. In: Trepińska, J., Olecki, Z. (Eds.), Klimatyczne aspekty środowiska
geograﬁcznego. Inst. Geogr. i Gosp. Przestrz. UJ, Kraków, pp. 139–148.
Żmudzka, E., 2002. O zmienności opadów atmosferycznych na obszarze Polski nizinnej w
drugiej połowie XX wieku. Wiadomości IMGW. 25 (46)(4), pp. 23–38.
Żmudzka, E., 2007. Zmienność zachmurzenia nad Polską i jej uwarunkowania cyrkulacyjne(1951-2000). Wydawnictwa Uniwersytetu Warszawskiego, Warszawa.
Żmudzka, E., 2009. Contemporary changes of climate of Poland [in Polish]. Acta Agrophys.
13 (2), 555–568.
Żurek, A., Mościcki, W., 2017. Application of penetrometer-based electrical resistivity proﬁling to study the unsaturated zone in a lysimeter experiment. Prace geograﬁczne 151,
121–132.

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, the Netherlands.
van Genuchten, M.T., 1980. A closed form equation for predicting the hydraulic conductivityof unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.
van Roosmalen, L., Christensen, B.S.B., Sonnenborg, T.O., 2007. Regional differences in climate change impacts on groundwater and stream discharge in Denmark. Vadose Zone
J. 6 (3), 554–571.
Von Hoyningen-Hüne, J., 1983. Die Interception des Niederschlags in landwirtschaftlichen
Beständen. Schriftenreihe des DVWK. 57, pp. 1–53.
Walker, G.R., 1998. Using soil water tracers to estimate recharge. In: Walker, G.R., Zhang, L.
(Eds.), The Basics of Recharge and Discharge. CSIRO Publication.
Wegehenkel, M., Kersebaum, K.C., 2009. An assessment of the impact of climate change on
evapotranspiration, groundwater recharge, and low-ﬂow conditions in a mesoscale catchment in Northeast Germany. J. Plant Nutr. Soil Sci. 172 (6), 737–744.
Wendland, F., Kunkel, R., Müller, U., Tetzlaff, B., 2008. Area – differentiated modelling ofgroundwater recharge rate for determining quantitative status of groundwater, casestudy: Federal State of Lower Saxony, Germany. In: Mikulič, Z., Andjelov, M. (Eds.),
Groundwater Modelling. Proceedings of Invited Lecturers of Symposium on Groundwater
Flow and Transport Modelling, Ljubljana, Slovenia, 28-31 January 2008, pp. 73–79.
Wibig, J., Glowicki, B., 2002. Trends of maximum and minimum temperature in Poland. Clim.
Res. 20, 123–133.
Wilcox, B.P., 2002. Shrub control and streamﬂow on rangelands: a process based viewpoint.
J. Range Manag. 55, 318–326.
Wine, M.L., Hendrickx, J.M., Cadol, D., Zou, C.B., Ochsner, T.E., 2015. Deep drainage sensitivity to climate, edaphic factors, and woody encroachment, Oklahoma,USA. Hydrol. Process. 29, 3779–3789.
Woldeamlak, S., Batelaan, O., De Smedt, F., 2007. Effects of climate change on the groundwater system in the Grote-Nete catchment,Belgium. Hydrogeol. J. 15, 891–901.
Wu, J., Zhang, R., Yang, J., 1996. Analysis of rainfall-recharge relationships. J. Hydrol. 177(1–2), 143–160.

15