Preliminary estimation of groundwater recharge on Brda river outwash plain

Estimation of groundwater recharge is one of the most challenging subjects in hydrogeology. It is a critical factor influencing the pollution migration, assessment of aquifer vulnerability to contamination, small-scale groundwater budget calculation, modeling of nutrient cycling and detailed flow path calculations. In Poland an infiltration rate method is widely used, which depends on a system of rate coefficients referring to the type of soil in the vadose zone and shows which part of the precipitation actually reaches the water table. The paper presents results of numerical simulations of vertical flow in unsaturated zone of an experimental site located on Brda outwash plain. Two simulations for varying vegetative cover (pine forest and grass) were performed. The results were compared with five infiltration rates classifications.


Introduction
Estimating the actual amount of recharge is critical in water management for the sake of protecting the quantity and quality of resources.Growing humans population and economic development increase water demand, of which main source is the groundwater [1].Lack of fresh water is a critical issue in many parts of the world, especially in regions where climate is arid or semiarid.Sustainable resources management, implemented in order to achieve rational and efficient consuming of water, makes groundwater recharge investigations essential to quantify safe yield of aquifer in hydrological balance of catchment.
Estimation of groundwater recharge is critical to contamination transport because it determines directions of groundwater flow in shallow aquifers [2].Therefore evaluation of infiltration rate is necessary for analysis of pollution migration, assessment of aquifer vulnerability to contamination, small-scale groundwater budget calculation and modeling, nutrient cycling, detailed flow path calculations [3].It is also important in waste disposal in order to location landfills outside of groundwater recharge zones to avoid contamination of water and soil with leakage [2].
Identifying recharge zones of wells, springs and hydrogeological bodies such as Major Groundwater Basins is required by Polish Environmental Law in order to protect water resources.There are many methods of assuming protection areas' sizes but most of them is based on vertical seepage time which is difficult to estimate [4].Correctness and accuracy of this assumption determine the cost of groundwater protection activities.
The main difficult in groundwater recharge estimation is the significant spatial and temporal variability of infiltration process.For groundwater modelers, planners and environmental policy makers recharge estimation is necessary, although in many cases they assign a single infiltration rate value for an entire watershed.This approach may result in serious errors in cases where detailed or small-scale flow paths have significant influence on the results.The amount of recharge is determined by many processes and phenomena which can be categorized into three groups: climatic, geomorphological (including topography, vegetation and soil) and geological factors [5].Such a number of factors make recharge rate one of the most complex hydrologic parameter to quantify [3].
The main goal of the present study was a preliminary estimation of groundwater recharge on an outwash plain of Brda river with a numerical model of vadose zone based on a representative soil profile.The results were compared with five widely used classifications of infiltration rates.The influence of vegetative cover (pine forest or grass cover) on the annual amount of recharge was also investigated.

Infiltration rate method
A number of techniques can be applied to estimate groundwater recharge [1,2].In humid regions, such as the European Plain, where precipitation significantly exceeds evapotranspiration, the main factor determining groundwater recharge is a subsurface geology, in particular ability of aquifer to store the water [2].Therefore in Poland an infiltration rate method is widely used [6].In this method recharge depend on system of rates which refer to the type of soil in vadose zone and represent the part of precipitation actually reaching the watertable (eg.[7]).
The infiltration rate method is the most common way of quantifying groundwater recharge in Poland.The coefficients of infiltration rate classifications reflect an average amount of water seeping to aquifer through different type of soil in the span of many years.Several classifications have been proposed, eg.Pazdro & Kozerski [7], Dyck & Chardabellas [8], Schneider & Züschang (cf.Załuski [9]), Wright et.al. [10] and Williams et al. [11].
Since classifications of Pazdro & Kozerski [7], Schneider & Züschang (cf.Załuski [9]) are generally applied in lowlands, classifications for uplands [12] and mountainous areas [13] have been introduced.Wright et.al. [10] and Williams et al. [11] classifications generate rather high infiltration rates as a result of specific climate of Ireland where they were implemented [6].Although this coefficients reflect the influence of soil type diversity quite well, none of them include impact of other factors aside from geology.During the works at Groundwater Vulnerability Map of Poland [14] a new methodology was introduced which considers main factors determining the infiltration process.The recharge rate estimation is based on equation: The method is recommended to calculations of well protection zones [15].

Numerical modeling
It is generally considered that, among models used to reflect water motion in unsaturated zone, Richards equation (RE) offers the most complete description of vertical flow in vadose zone [1] with all microhydrologic processes such as capillary rise, dewatering of soil, retention of water in soil, root water uptake, evaporation and percolation.RE couples Darcy law with mass conservation equation as follows: where: . Hydraulic conductivity is usually described as a function of negative pressure head [K(h)] or as a function of water content [K()].RE has a form of nonlinear partial differential equation because of coefficient K() and K(h) which are strongly nonlinear functions [16] and it is recommended to identify them empirically way.
In our study we solved the Richards equation using HYDRUS-1D computer code [17] which simulates the transient movement of water, solute, and heat in variably-saturated soil profiles using finite element discretization in space and fully implicit discretization in time.The software has been validated in a large number of applications, including several related to recharge (e.g.[18,19]).

Site description
The area of investigation is located in the Pomerania Region (Northern Poland), in Cekcyn, 10 km south-east from the town of Tuchola.It is a part of a young-glacial Brda river outwash plain, which combined with the adjacent outwash plain of Wda river, is perceived as the largest hydrosystem of this type in Poland [20].The upper, unconfined aquifer is located in fluvio-glacial sands of the outwash plain.The depth to water table is about 6-9 m and the aquifer thickness is about 5-7 m.Below the upper aquifer is a layer of glacial till (moraine), about 5-9 m thick, separating it from the lower (confined) aquifer.The piezometric head in the lower aquifer is slightly (0.6-1 m) below the water table in the unconfined aquifer.Therefore, the upper aquifer is recharged only by infiltration, without any lateral inflow or inflow through the clay layer.Four measurement profiles (PR1-PR4) was installed at the site, in order to account for variability in vegetation (grass cover vs. coniferous (pine) forest).In each profile soil water content at 4 different depths (20 cm, 50 cm, 150 cm, 450 cm) and the level of groundwater table has been measured continuously since April 2017, with remote data logger.Hydraulic conductivity for each layer, investigated with field permeameter, varies from 0.6 m/d (sandy loam) to 30.5 m/d (medium sand).Results of the measurements are shown on fig. 2.

Model setup
Two simulations were performed basing on one dimensional model of representative vadose zone profile PR4.The first scenario referred to a grassy area whereas the second one represents pine woodland.Some geological and hydrogeological parameters of the model, such as unsaturated zone thickness, vertical variety of hydraulic conductivity and water content in soil were measured during field investigation, weather data was obtained from nearby meteorological station in Chojnice and the remaining data, e.g.parameters of root water uptake were assigned according to literature [21][22][23].Simulations were carried out for the years 2003-2016 with daily time step.
The model has 7 meters of thickness and includes five different soil materials (fig 2).Vertically it is divided with 71 nodes, placed with larger intervals in the lower, deeper part of the model and smaller distances in the shallow subsurface, in order to better reflect the variability of flow conditions at the soil atmosphere interface.The upper boundary condition represents varying in time precipitation (with surface runoff) while the lower boundary condition was assigned as a constant pressure head.
The retention function and relative conductivity function was described with the Brooks-Corey model [24], with representative parameters for sand and sandy loam, taken from [21] (Table 1).

Fig. 2. Parameters of numerical model.
Weather data required to simulation (annual temperature, air humidity observations and amount of precipitation) was obtained from Chojnice meteorological station.Potential evapotranspiration was calculated with Grabarczyk formula [25].The obtained values were distributed proportionally through the root zone as a linear function which decline with depth.Pine and grass actual root water uptake was estimated with Faddes' model [26] based on parameters presented in references [22,23].The thickness of the plant root zone was assigned according to Meyer et al. [27] as 2 m for pine forest and 0.5 m for grass cover.

Results and discussion
Results show that there is a significant difference in groundwater recharge between both case scenarios (fig.3).During 14 years of simulation aquifer received average annual 302 mm/y (49% of rainfalls) on forested area and average annual 332 mm/y (53% of precipitation) on grassland.The results were confronted with five classifications of infiltration rates (Table 2).In general the methods dedicated to estimations of groundwater recharge on lowlands range from 0.27 to 0.38classifications of Dyck & Chardabellas [8] and Duda et al. [14] give similar values depending on vegetative cover while Pazdro & Kozerski [7] system provides average coefficient of recharge.Comparing to this values, results of numerical model seems to be excessive.Rather high infiltration rates of Wright [10] and Williams et al. [11] classifications are perceived as a result of developing this method in Ireland where climate is different than in Poland [6].Overestimated infiltration rates obtained from the simulations can be caused by relatively low values of evapotranspiration.Some studies shows that Grabarczyk formula provide underestimated results compared with different models (eg.Penmana-Monteith's, Hargreaves') [28].The results of numerical simulations correlated with classifications of Dyck & Chardabellas [8] and Duda et al. [14] showed that vegetative cover has an effect on groundwater recharge.The actual root water uptake in the numerical model ranges from 291 mm/y in grassland to 325 mm/y for pine forest.The differences between the amount of percolating water for pine woodland and grassland reach 4% of precipitation based on simulation results, 8% according to method of Dyck & Chardabellas [8] and 9% in Duda et al. [14] classifications.Plants can absorb significant part of precipitation for their life processes.The root zone plays the role of a buffer which reduces the amount of infiltrating water and allows only some part of it to drain towards deeper groundwater table [29].However, this is only a preliminary evaluation of recharge rate, which needs to be verified in detail with future field observations, a number of studies show that plant cover influence groundwater recharge in a significant way [23] not only in point-scale estimations but also in regional assessment [3].

Table 1 .
[24]meters of the Brooks-Corey[24]hydraulic model for the soils used in simulations.