SoftwareX 26 (2024) 101680

Contents lists available at ScienceDirect

SoftwareX
journal homepage: www.elsevier.com/locate/softx

Original software publication

HMSE: A tool for coupling MODFLOW and HYDRUS-1D computer programs
Mateusz Pawlowicz a , Bartosz Balis a ,∗, Adam Szymkiewicz b , Jirka Šimůnek c ,Anna Gumuła-Kawęcka b , Beata Jaworska-Szulc b
a AGH University of Krakow, Krakow, Polandb Gdansk University of Technology, Gdansk, Polandc University of California-Riverside, CA, USA

ARTICLEINFO

Keywords:
Groundwater
HYDRUS-1D
MODFLOW
Model coupling
Kubernetes

ABSTRACTA new software HMSE has been developed to facilitate external coupling between two well-known programs for
subsurface flow modeling: MODFLOW-2005 (saturated zone flow) and HYDRUS-1D (unsaturated zone flow).
Two coupling schemes have been implemented. In the first case the groundwater recharge flux is calculated by
HYDRUS-1D assuming a fixed water table position and then passed to MODFLOW input files. In the second casethe water table position in HYDRUS-1D is updated periodically using the solution from MODFLOW. HMSE canbe deployed in 3 modes: local, Docker and Kubernetes cluster. A web-based interface is provided to configure
and run the simulation in all three cases. The software is applied to simulate groundwater table fluctuations
observed in a shallow aquifer during three years.

Code metadata
Current code version
Permanent link to code/repository used for this code version
Code Ocean compute capsule
Legal Code License
Code versioning system used
Software code languages, tools, and services used
Compilation requirements, operating environments & dependencies
If available Link to developer documentation/manual
Support email for questions

1. Motivation and significanceNumerical models are increasingly used for groundwater management. They allow us to understand better the past, present, and future
conditions of groundwater systems, find optimal management scenarios, implement protection measures for wellheads and aquifers, andevaluate the impact of climate and land use change on groundwater
resources. Water in the subsurface occurs in two distinct zones: theunsaturated zone (or vadose zone) between the ground surface andgroundwater table and the saturated zone (groundwater zone) below

1.0.0
https://github.com/ElsevierSoftwareX/SOFTX-D-23-00705
MIT License
Git
Python, Javascript, Kubernetes, Airflow, MinIO
Python >= 3.8, Airflow >= 2.3.0
https://github.com/WaterlinePL/HMSE-frontend#readme
adam.szymkiewicz@pg.edu.pl
balis@agh.edu.pl
mat.pawlowicz1@gmail.com

the groundwater table. In the unsaturated zone the pores of soils androcks are filled partly with water and partly with air, while in thesaturated zone the pores are completely filled with water. From thecomputational point of view, the simulation of flow in the vadosezone is much more challenging than in the saturated zone, because
the governing equations are highly nonlinear and require a significant
number of parameters which are often difficult to obtain.
Groundwater is extracted from aquifers, highly permeable geological formations in the saturated zone. However, the quantity and qualityof water in aquifers critically depends on the processes occurring in theunsaturated zone. Groundwater is replenished mainly by water from∗ Corresponding author.

E-mail address: balis@agh.edu.pl (Bartosz Balis).
https://doi.org/10.1016/j.softx.2024.101680
Received 17 October 2023; Received in revised form 27 February 2024; Accepted 4 March 2024
2352-7110/© 2024 Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

SoftwareX 26 (2024) 101680

M. Pawlowicz et al.

precipitation which seeps through the unsaturated zone and ultimately
arrives at the groundwater table as recharge. The recharge amount isstrongly affected by climatic conditions, the permeability of soils in
the unsaturated zone, the type of plants and the depth of groundwatertable. Recharge determines the amount of groundwater that can besustainably extracted from aquifers and understanding its spatial andtemporal variability is crucial for efficient groundwater management.
If the groundwater table is shallow, water can be taken up by capillaryaction and lost due to evapotranspiration (combined action of evaporation and uptake by plants). Thus, accurate representation of vadosezone processes can significantly improve groundwater modeling.
A number of integrated hydrological models have been developedthat consider both unsaturated and saturated zones as a 3D continuum(e.g. PAR-FLOW [1], CATHY [2], Hydro-Geo-Sphere [3], FEFLOW [4],HYDRUS (2D/3D) [5]). However, their use in routine groundwater
management applications is still limited, because of their high demandfor input data, computer resources, and user expertise. Groundwatermodels used in management practice are commonly developed using
the MODFLOW family of programs [6,7], which simulate 3D waterflow in the saturated zone. MODFLOW can be extended by couplingwith 1D vadose zone models, describing water flow in representative
soil profiles associated with selected parts of the groundwater model
domain. This is justified as the flow in the vadose zone occurs mostly
in the vertical direction. This approach has several advantages: (i) it
is more efficient computationally than a full 3D saturated/unsaturated
model (1D simulations of the vadose zone can be easily parallelized);(ii) it allows for representing a variety of vadose zone processes relevant to groundwater flow in realistic, physically-based manner (1Dmodels were shown to be useful in representing spatial and temporal variability of recharge and predicting the effects of land use
and climate change (e.g. [8])); (iii) vadose zone modeling can beapplied only to selected parts of the groundwater model domain. The
UZF package [9] included in MODFLOW simulates 1D gravity-driven
flow in the vadose zone using the kinematic wave approach. While
it is very efficient computationally, the kinematic wave model neglects capillary-driven flow, which can be important in some situations(e.g. upward flow caused by root water uptake). The Richards equation
provides a more comprehensive vadose zone flow model, which considers both gravity and capillary forces. State-of-the art simulators basedon the Richards equation are available, including HYDRUS-1D [10] andSWAP [11] and have been coupled to MODFLOW programs [12–14].
The coupling of HYDRUS-1D and MODFLOW was initially done for
MODFLOW-2000 and later applied to MODFLOW-2005 [15–17]. It wasalso extended to solute transport modeling with MT3DMS [18]. However, the previous versions were coupled on the code level, i.e. parts ofthe computer code of HYDRUS-1D were merged into MODFLOW (bothare written in Fortran) and a new standalone computer application was
created. One disadvantage of such an approach is that it cannot be usedwith standard input files of HYDRUS-1D and MODFLOW and creating
application-specific input files is awkward. Moreover, any modification
or update of standard MODFLOW or HYDRUS-1D codes must be ported
manually into the merged code.
In this work we developed a new coupling approach for HYDRUS-1D
and MODFLOW-2005: HYDRUS-MODFLOW Synergy Engine (HMSE).
We use external coupling via an independent application, which allows
the user to prepare MODFLOW and HYDRUS-1D models independently,
using the graphical user interfaces available for both programs. The
application modifies input files to pass information between HYDRUS1D and MODFLOW, but each of these two programs is called externallyas an executable file.

existing open-source Python libraries for handling MODFLOW models(FloPy, [19]) and HYDRUS-1D (PHydrus, https://github.com/phydrus/
phydrus). The following sections briefly describe the details of thesoftware. More detailed benchmark of the software on a simple model
can be found in the Supplementary material.
2.1. Software architecture
The logical architecture of the software is shown in Fig. 1 . The
yellow and purple blocks describe system elements present in all the deployments. The dashed boxes indicate components which differentiate
among application versions. The orange blocks describe components
used in desktop deployment, the blue ones are used by Docker version
and the gray ones are used by Kubernetes deployment. Users interact
with a Web User Interface, a web application powered by a Python
WSGI HTTP Server. The interface allows for setting up projects whereusers can define recharge zones, upload HYDRUS-1D and MODFLOW
models, set various configuration parameters, run the simulations, anddownload the results.
The project’s configuration is handled by the Project configuration
module, which facilitates the setup of hydrological simulations. The
projects’ storage depends on the type of deployment. In a desktop
deployment, this module interfaces with the local file system. In Docker
deployment, it uses container volumes mapped to the host file system.
Finally, in a Kubernetes cluster deployment, it communicates with an
S3-compatible object storage (e.g., MinIO) through its API.
The Scheduling component supervises hydrological simulations, orchestrates actions and invokes simulations using third-party programs(HYDRUS-1D and MODFLOW-2005). Both the desktop deployment andthe Docker version create Python supervisor threads for this purpose.
By contrast, the Kubernetes cluster version leverages the Airflow scheduler for running the simulation pipeline as a set of Kubernetes Jobs ina cluster of machines.
Finally, the Workflow environment is where all the simulations
occur. In the desktop version, the simulations are performed on copies
of project files within the local file system by creating OS processes
to run simulations. The Docker deployment operates similarly, butusing a shared volume and the Docker-in-Docker technique to run thesimulation steps as containers. On the other hand, the Kubernetes deployment runs the HYDRUS models in parallel using Pods in a dedicatednamespace. Consequently, a distributed file system is needed to sharedata between the simulation steps. This can be achieved, for example,
through an NFS Persistent Volume (PV) provided by the cluster. The
simulation pipeline first retrieves the simulation data files from an
S3 storage and places them in a working directory on the shared filesystem. Then the simulation steps are executed as Kubernetes Jobs
(Pods). Finally, the results are uploaded back to S3.
2.2. Software functionalities
The coupling between MODFLOW and HYDRUS-1D models is implemented in HMSE along the lines described in [12,13]. A conceptual
scheme is shown in Fig. 2. The area of the MODFLOW model is divided
into several recharge zones, and HYDRUS-1D models are assigned toall or selected zones. Each HYDRUS-1D model represents a typicalsoil profile, land cover, depth to groundwater and weather conditionsin the considered zone. MODFLOW model then represents transientgroundwater flow, the simulation of which is divided into several
user-defined stress periods, representing variations in boundary conditions, groundwater extractions etc. The main information obtained
from HYDRUS-1D and passed to MODFLOW is groundwater recharge
in each stress period. HMSE offers two coupling modes. In the simple
coupling mode the groundwater table position is assumed constant intime and corresponding to the bottom of each HYDRUS-1D profile.
First, simulations for all HYDRUS profiles are performed for the entire
considered time period. Then, the average recharge rate corresponding

2. Software description
The HMSE comes in 3 deployments: desktop application (standalone
webserver), docker container (standalone webserver) and Kubernetes
cluster (microservice architecture). A part of HMSE code is based on
2

SoftwareX 26 (2024) 101680

M. Pawlowicz et al.

Fig. 1. Logical architecture of HMSE. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)then the recharge data is passed to the RCH file and a MODFLOW
simulation is run for a single stress period. Based on the MODFLOW
results, the average depth to the water table is calculated for each
recharge zone and the corresponding HYDRUS profile is updated. The
algorithm for updating water pressure distribution in the soil profile
allows for avoiding spurious oscillations in the recharge flux due todiscontinuities in the water table position between subsequent stressperiods. It is described in detail in [16].
HMSE runs using the default web browser on a given machine. In
the Configure tab the user has to specify paths to MODFLOW-2005
and HYDRUS-1D executables, which must be independently installedon the computer. Work in HMSE is organized into projects. Each projectcontains a single MODFLOW-2005 model (the user must upload allmodel files as a single *.zip archive). A number of models developedin HYDRUS-1D for various parts of the MODFLOW model domain
can be added by the user (each HYDRUS model is uploaded as asingle *.zip archive containing all input text files). Additionally, atleast one file with weather data must be added. The weather file is in
CSV format and contains daily precipitation, maximum and minimum
temperature, relative humidity, wind speed, solar radiation, as well as
the latitude, longitude and elevation of the weather station. The user
assigns a weather file to each HYDRUS profile and a HYDRUS model
to one or more recharge zones in the MODFLOW model. It should be
noted that weather data are typically provided in HYDRUS-1D input
files atmosph.in and meteo.in. In HMSE these files are overwrittenwith data from the weather files provided by the user as a separate
input. The main purpose of such solution was to facilitate the use ofdifferent weather data sets and different simulation periods with thesame HYDRUS-1D profiles. Weather data can be changed without theneed for the user to update manually the HYDRUS profiles, which isuseful especially in long-term simulations. Potential evapotranspiration
is calculated internally by HYDRUS-1D, based on the weather data andother parameters specified by the user in the definition of soil profile(such as leaf area index LAI or crop height). In the future we plan tosupport more input options for weather data, including user-defined
potential evapotranspiration.
HMSE offers two ways to define recharge zones in the MODFLOW
model. The first method is based on the analysis of the existing RCH
file in the model provided by the user. Various recharge zones are
typically defined in existing groundwater models, based on geological,hydrological, and other factors. The cells in a single zone are assigned
the same value of a recharge flux in a given stress period in the
RCH file. HMSE parses the RCH file of the original MODFLOW modeland divides the model domain into several recharge zones. In theFig. 2. Conceptual scheme of coupling between HYDRUS-1D and MODFLOW
implemented in HMSE.

to each zone in each stress period is calculated by HMSE and written
to the MODFLOW RCH file, which provides recharge rates for the
groundwater flow simulation. In the second stage, the MODFLOW
simulation is carried out.
The second coupling mode includes feedback from MODFLOW toHYDRUS-1D, to update the position of the water table in HYDRUS1D profiles after each stress period in the MODFLOW simulation. To
achieve this, the MODFLOW model is split into a number of models
corresponding to subsequent stress periods of the main simulation. For
each stress period the HYDRUS-1D simulations are performed first,
3

SoftwareX 26 (2024) 101680

M. Pawlowicz et al.

Fig. 3. HMSE interface screenshot showing the assignment of recharge zones.

second method the recharge zones correspond to the zones definedwith the ZONEBUDGET [20] tool, which is often used as a postprocessor to MODFLOW. In this case HMSE reads the text file defining
ZONEBUDGET zones, which must be provided by the user together
with other MODFLOW files in the *.zip archive. The preprocessing
functionalities of HMSE are limited to handling the recharge zones
and weather data. The MODFLOW and HYDRUS-1D models must be
prepared independently by the user using other available tools, such as
GMS,1 ModelMuse [21], or FloPy [19] (for MODFLOW-2005) and the
HYDRUS-1D GUI2 or PHydrus library3 (for HYDRUS-1D). For running
the HYDRUS models under HMSE we recommend to use the HYDRUS1D versions available in the GitHub repository for Windows and Linuxsystems, which are based on the open-source version of HYDRUS-1D,4
with some minor changes to improve stability and ensure compatibility
with the most recent version of HYDRUS-1D (v. 4.17). Currently HMSE
is only compatible with MODFLOW-2005 models, however work is inprogress to extend coupling to MODFLOW 6.
Once the recharge zones are defined, the uploaded HYDRUS-1Dmodels can be assigned to one or more zones, Fig. 3. The user also
defines the starting date of the simulation period, which is used to
extract weather data from the external file. It is also possible to define
a warm-up period. In the warm-up period HYDRUS-1D models are runwith the corresponding weather data. The water pressure distribution
at the end of the warm-up period is then used as a realistic initialcondition for HYDRUS simulations in the time period of interest.

isotropic hydraulic conductivity 𝑘. From top to bottom the layers were:(1) outwash sand and gravel (upper aquifer, 𝑘 = 25 m/d), (2) glacial
till (𝑘 = 0.05 m/d), (3) outwash sand (lower aquifer, 𝑘 = 10 m/d), (4)silt and glacial till (𝑘 = 0.01 m/d). Water flow between the surface
water (lakes) and the aquifers (layers (1) and (3)) was represented
with a third-type boundary condition (General Head Boundary). In theaquitard layers (2) and (4) a no-flow condition was imposed at lakeboundaries. The northern boundary was perpendicular to groundwaterhead contours, so it was assigned a no-flow (second-type) condition.
Two recharge zones with different plant cover (pine forest andgrassland) were distinguished (Fig. 4). The two corresponding soil
profiles were taken from [22] and are denoted here as P1 and P3,
consistent with that study. Each profile had two sandy layers (Sand I
and II) separated by a layer of sandy loam. The hydraulic characteristics
of each soil material were described with the van Genuchten–Mualem
model [24]. An atmospheric boundary condition with no water ponding
option (instantaneous runoff) was specified on the soil surface basedon daily weather data measured on-site. Potential evapotranspiration(PET) was calculated using the Penman–Monteith equation [25] witha leaf area index (LAI) equal to 3.0 in the pine forest and 2.0 in thegrassland. The vegetation height was set to 30 cm (it is not possible
to specify realistic values for forests with the current implementationof PET calculation, nevertheless, the resulting PET estimates werereasonable for both land covers). Root water uptake was estimatedusing the Feddes macroscopic model [26] with a root zone depth of
1.5 m for the forest and 0.5 m for the grass cover. The reader is referredto [22] for a detailed soil characterization and parameters. Simulations
were performed for a 3-year period (15 April 2017–14 April 2020),
with a 3-year warm-up simulation carried out to determine initial flow
conditions. HYDRUS-1D used automatically adjusted time step with amaximum of 1 day. MODFLOW simulation was divided into 36 monthly
stress periods.
Simulations were carried out using the two coupling methods described above. In the first scenario (one-way coupling) profiles were
6.60 m (forest), and 7.05 m (grassland) deep and the bottom boundary
condition was assumed to be a constant pressure head equal to zero(water table level). In the second scenario (two-way coupling), theprofiles were extended to 8 m, and the water pressure head at thebottom was updated according to the results obtained from MODFLOW.
The recharge fluxes obtained using the two coupling methods areshown in Fig. 5. The one-way coupling results in a smoother evolution

3. Illustrative example
HMSE (desktop version) has been applied to model a shallow aquifer
on an outwash plain in northern Poland, as studied in [22]. The aquifer
has horizontal extent of about 2 km2 and is enclosed by four lakes
in the west, the south and the east and a partially wet depression
in the north (Fig. 4). The MODFLOW model consisted of gridblocks(cells) with horizontal dimensions of 10 by 10 m arranged in four
layers (74 770 cells in total). Each layer was assigned a uniform and1
https://www.aquaveo.com/software/gms-groundwater-modeling-systemintroduction
2
https://www.pc-progress.com/en/Default.aspx?hydrus-1d3
https://github.com/phydrus/phydrus)
4https://www.pc-progress.com/en/Default.aspx?H1D-description#k8

4

SoftwareX 26 (2024) 101680

M. Pawlowicz et al.

Fig. 4. Model area and characteristic soil profiles.
Source: Modified from [23].

4. Impactof recharge in time, while two-way coupling leads to more fluctuations due to the changes in the position of the water table between
MODFLOW stress periods. Negative recharge values appear in two-way
coupling, indicating upward water movement due to capillary uptake.

The current version of HMSE allows for extending existing groundwater models developed with MODFLOW-2005 by coupling them with
HYDRUS-1D models simulating hydrological processes in representative vadose zone profiles. Both HYDRUS-1D and MODFLOW-2005 aremature, state-of-the-art computer programs validated on thousandsof case studies and supported by large user communities. Our application brings synergy effects of more comprehensive hydrologicalmodeling based on the combined use of these two programs. Including detailed description of vadose zone processes can be especiallyhelpful in: (i) improved understanding of time and space variability ofgroundwater recharge, (ii) establishing protection zones for groundwater intakes and aquifers, (iii) evaluation of the impact of climate andland use change on groundwater resources, (iv) assessing vulnerability of aquifers to pollution originating from the ground surface, (v)evaluation of managed aquifer recharge (MAR) schemes, which areoften based on artificial infiltration of water from the surface via thevadose zone to aquifers. Future work will focus on the implementationof coupling for solute transport and integration with MODFLOW 6

The simulated water table elevation was compared to field measurements carried out in P1 and P3 [22] (Fig. 6). The trends andmagnitude of water table movement in simulations show reasonablesimilarity to the observations, even though the agreement is less thanperfect. The differences are larger for the grassland than for the forest.
Also, the one-way coupling seems to give the results more alignedwith observations than the two-way coupling, especially during the dry
period, when the simulation overestimates the magnitude and speed of
the water table drawdown. This might be attributed to several reasons,such as variability in aquifer parameters, extrapolating data from single
soil profiles to larger areas and comparing point measurements of thegroundwater table to area averages obtained from MODFLOW. The
fit could possibly be improved by recalibrating the coupled model.
However, such a task is beyond the scope of our study.
5

SoftwareX 26 (2024) 101680

M. Pawlowicz et al.

Fig. 5. Groundwater recharge flux calculated by HYDRUS.

Fig. 6. Observed and simulated groundwater table elevation.

5. Conclusions

Acknowledgments

A new application has been developed that enables to couple
MODFLOW-2005 groundwater models with HYDRUS-1D vadose zone
models in a user-friendly way. The application allows two coupling
modes and two methods to delineate recharge zones in MODFLOW
models. Preliminary validation has been performed on a 3-year seriesof observations in a shallow sandy aquifer and showed satisfactory
performance of the new computer tool.

Development of the HMSE software was funded by the project WATERLINE, which is supported by the National Science Centre, Poland(2020/02/Y/ST10/00065), under the CHIST-ERA IV programme, which
has received funding from the EU Horizon 2020 Research and Innovation Programme, under Grant Agreement no 857925.

CRediT authorship contribution statement

Supplementary material related to this article can be found online
at https://doi.org/10.1016/j.softx.2024.101680.

Appendix A. Supplementary data

Mateusz Pawlowicz: Writing – review & editing, Writing – originaldraft, Software. Bartosz Balis: Writing – review & editing, Writing –original draft, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Conceptualization. Adam Szymkiewicz:
Writing – review & editing, Writing – original draft, Supervision, Resources, Methodology, Investigation, Conceptualization. Jirka Šimůnek:Software, Conceptualization. Anna Gumuła-Kawęcka: Writing – review & editing, Writing – original draft, Visualization, Validation,
Resources, Investigation, Conceptualization. Beata Jaworska-Szulc:
Writing – review & editing, Writing – original draft, Validation, Resources, Conceptualization.

References[1] Maxwell RM, Kollet SJ, Smith SG, Woodward CS, Falgout RD, Ferguson IM, et al.
Parflow user’s manual. Tech. rep., International Ground Water Modeling Center
Report GWMI; 2009.
[2] Camporese M, Paniconi C, Putti M, Orlandini S. Surface-subsurface flow modeling with path-based runoff routing, boundary condition-based coupling, andassimilation of multisource observation data. Water Resour Res 2010;46(2).
[3] Therrien R, McLaren R, Sudicky E, Panday S. Hydrogeosphere: a threedimensional numerical model describing fully-integrated subsurface and surface
flow and solute transport. Tech. rep., Waterloo, ON: Groundwater Simulations
Group, University of Waterloo; 2010.
[4] Diersch H-JG. FEFLOW: finite element modeling of flow, mass and heat transportin porous and fractured media. Springer Science & Business Media; 2013.
[5] Šejna M, Šimůnek J, van Genuchten MT. The HYDRUS software package forsimulating the one-, two- and three-dimensional movement of water, heat andmultiple solutes in variably-saturated porous media. User manual ver. 5.x. Tech.
rep., Prague, Czech Republic: PC-Progress; 2022.
[6] Harbaugh AW, Banta ER, Hill MC, McDonald MG. MODFLOW-2000, the US
geological survey modular ground-water model-user guide to modularization
concepts and the ground-water flow process. VA, USA: USGS Reston; 2000.
[7] Harbaugh AW. MODFLOW-2005, the US geological survey modular ground-watermodel: the ground-water flow process, vol. 6, VA, USA: USGS Reston; 2005.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared toinfluence the work reported in this paper.
Data availability
Data will be made available on request.
6

SoftwareX 26 (2024) 101680

M. Pawlowicz et al.

[17] Szymkiewicz A, Gumuła-Kawęcka A, Šimůnek J, Leterme B, Beegum S, JaworskaSzulc B, et al. Simulations of freshwater lens recharge and salt/freshwaterinterfaces using the HYDRUS and SWI2 packages for MODFLOW. J Hydrol
Hydromech 2018;66(2):246–56.
[18] Beegum S, Šimůnek J, Szymkiewicz A, Sudheer K, Nambi IM. Implementation ofsolute transport in the vadose zone into the HYDRUS package for MODFLOW.
Groundwater 2019;57(3):392–408.
[19] Bakker M, Post V, Langevin CD, Hughes JD, White JT, Starn J, et al.
Scripting MODFLOW model development using python and flopy. Groundwater
2016;54(5):733–9.
[20] Harbaugh AW. A computer program for calculating subregional water budgets
using results from the US geological survey modular three-dimensional finitedifference ground-water flow model, vol. 90, VA, USA: US Geological Survey
Reston; 1990.
[21] Winston RB. Modelmuse: a graphical user interface for MODFLOW-2005 andphast. Tech. rep., VA, USA: US Geological Survey Reston; 2009.
[22] Gumuła-Kawęcka A, Jaworska-Szulc B, Szymkiewicz A, Gorczewska-Langner W,
Pruszkowska-Caceres M, Angulo-Jaramillo R, et al. Estimation of groundwater
recharge in a shallow sandy aquifer using unsaturated zone modeling and watertable fluctuation method. J Hydrol 2022;605:127283.
[23] Gumuł a Kawęcka A, Szymkiewicz A, Jaworska-Szulc B, Pruszkowska-Caceres M,
Gorczewska-Langner W. Preliminary estimation of groundwater recharge on
brda river outwash plain. In: 10th conference on interdisciplinary problems inenvironmental protection and engineering EKO-DOK 2018. EDP Sciences; 2018,
p. 00050.
[24] van Genuchten MT. A closed form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci Am J 1980;44(5):892–8.
[25] Allen R, Pereira L, Raes D, Smith M. Crop evapotranspiration-guidelines forcomputing crop water requirements-FAO irrigation and drainage paper 56. Tech.
rep., Rome: FAO; 1998.
[26] Feddes R, Kowalik P, Zaradny H. Simulation of field water use and crop yield.
New York, NY: John Wiley and Sons; 1978.

[8] Gumuła-Kawęcka A, Jaworska-Szulc B, Szymkiewicz A, Gorczewska-Langner W,
Angulo-Jaramillo R, Šimůnek J. Impact of climate change on groundwater
recharge in shallow young glacial aquifers in northern Poland. Sci Total Environ2023;877:162904.
[9] Niswonger RG, Prudic DE, Regan RS. Documentation of the unsaturated-zone
flow (uzf1) package for modeling unsaturated flow between the land surface
and the water table with MODFLOW-2005. Tech. rep., USGS; 2006.
[10] Šimůnek J, Šejna M, Saito H, SM, Van Genuchten MT. The HYDRUS-1D software
package for simulating the one-dimensional movement of water, heat, andmultiple solutes in variably-saturated media. Ver. 4.17. Tech. rep., Departmentof Environmental Sciences, University of California-Riverside; 2013.
[11] Kroes J, Van Dam J, Bartholomeus R, Groenendijk P, Heinen M, Hendriks R,
et al. Swap version 4: theory description and user manual. Tech. rep.,
Alterra-rapport-Wageningen University and Research Centre; 2017.
[12] Seo HS, Šimůnek J, Poeter EP. Documentation of the HYDRUS package for
MODFLOW-2000, the us geological survey modular ground-water model. Tech.
rep., IGWMC-International Ground Water Modeling Center; 2007.
[13] Twarakavi NKC, SŠimůnek J, Seo S. Evaluating interactions between groundwaterand vadose zone using the HYDRUS-based flow package for MODFLOW. Vadose
Zone J 2008;7(2):757–68.
[14] Xu X, Huang G, Zhan H, Qu Z, Huang Q. Integration of SWAP and MODFLOW2000 for modeling groundwater dynamics in shallow water table areas. J Hydrol
2012;412:170–81.
[15] Leterme B, Gedeon M, Jacques D. Groundwater recharge modeling of the Netecatchment (Belgium) using the HYDRUS-1D–MODFLOW package. In: Proceedingsof the 4th international conference HYDRUS software applications to subsurface
flow and contaminant transport problems. 2013, p. 235–43.
[16] Beegum S, Šimůnek J, Szymkiewicz A, Sudheer K, Nambi IM. Updating thecoupling algorithm between HYDRUS and MODFLOW in the HYDRUS package
for MODFLOW. Vadose Zone J 2018;17(1):1–8.

7