Dominant Influencing Factors of Groundwater Recharge Spatial Patterns in Ergene River Catchment, Turkey

Abstract

Groundwater is of great significance in sustaining life on planet earth. The reliable estimation of groundwater recharge is the key understanding the groundwater reservoir and forecasting its potential accessibility. The main objective of this study was to assess the groundwater recharge and its controlling factors at the Ergene river catchment. A grid-based water balance model was adopted to determine the spatially distributed long-term groundwater recharge and other water budget components, relying upon the hydro-climatic variables, land-use, soil, geology, and relief of the investigated area. The model calculations were performed for the hydrological reference horizon of 20 years at a spatial resolution of 100 × 100 m. The base flow index (BFI) separation concept was applied to split up the simulated total runoff into groundwater recharge and direct runoff. Subsequently, the statistical methods of Pearson product–moment correlation and principal component analysis (PCA) were combined for identifying the dominant catchment and meteorological factors influencing the recharge. The average groundwater recharge over the investigated area amounts to 95 mm/year. The model validation and statistical analysis indicate that the difference between simulated and observed total runoff and recharge values is generally under 20% and no significant inconsistency was observed. PCA indicated that recharge is controlled, in order of significance, by land-use, soil, and climate variables. The findings of this research highlight the key role of spatial variables in recharge determination. In addition, the generated outputs may contribute to groundwater resource management in the Ergene river catchment.

1. Introduction

Groundwater is the capital source of freshwater for nearly half of the earth’s population for irrigation and domestic water needs [1]. It is identified as a renewable water resource for supporting agricultural, industrial, environmental, and municipal domestic water demands [2]. The estimation of groundwater recharge is the key to understanding the groundwater reservoir and forecasting its potential accessibility and sustainability, even though other elements have to be taken into account; e.g., social, economic, and hydrogeological considerations [3]. Nevertheless, groundwater recharge is the least understood component, mainly due to its extensive spatial and temporal variation and that it requires much effort for direct measurement [4]. The studies about groundwater recharge in a catchment or on a regional scale are crucial for determining the quantitative status of groundwater resources.

Existing studies about groundwater in the Ergene catchment focused only on the groundwater quality [5,6], not on the quantity. The spatially distributed water budgets, including the recharge, have not been broadly ascertained and properly understood. Despite this, the groundwater is the central and vital source of municipal domestic water for 1,150,000 people and industrial water demand for 2500 leather and textile businesses; furthermore, agriculture covers 73% of the total land of the Ergene catchment [7]. Sustainable water resource management, where the influence of agricultural activity over a catchment area is high, aims to guarantee a durable, convenient, and long-standing water supply system to achieve the plants’ water needs, along with an increase in a number of industries and municipal drinking water needs, through mitigating harmful environmental effects. Fundamental water resource management is directly linked with the estimation of accessibility of water within a drainage area and implies the understanding of catchment groundwater hydrologic processes. In this regard, the spatial-temporal pattern of the groundwater recharge generated in this study will have a significant role in supporting the water resource management scheme of the Ergene catchment. Groundwater recharge refers to deep percolation of water from the land surface through the vadose zone and replenishing of the groundwater aquifer. Its quantitative assessment plays a central role in the awareness of large-scale hydrological processes and groundwater supply management and sustainability. To expand the approach of water resources management in the Ergene River catchment, the water budget components, e.g., groundwater recharge and direct runoff, have to be separated and calculated. Estimation of groundwater recharge rates are indispensable for sustainable groundwater management and proper planning in the Ergene catchment.

Various methods may be used to estimate the recharge [8], including direct measurement techniques, e.g., seepage meters and tracer methods. However, the field measurement techniques seem to be extremely hard [9], and are influenced by assessment errors, spatial irregularity patterns, and are operated at high cost [10]. The recorded baseflow from the river gauging stations and its analysis can be used to estimate the level of groundwater recharge [11,12] as a comparison method. Base flow is taken into account as a groundwater recession flow, which is the portion of groundwater flow to streams during rainless seasons. The main consideration in this approach, in order to determine the recharge levels, is that base flow is taken as equivalent to the total groundwater recharge of a basin and the groundwater discharge is roughly considered as the recharge [13]. The base flow analysis method does not consider the origin of the water; however, it considers that great rates of base flow are generated by flow coming from the groundwater reservoir. Thus, for a long period of time, base flow rates were generally taken as watershed groundwater recharge [14]. This method is broadly useful to validate groundwater models which are based on empirical water balance techniques [15]. However, this approach involves the complexity of directly comparing the recharge and base flow, while the latter approach estimates the substitution of groundwater discharge and thus of real recharge [13]. Computer-based models, which rely on a soil–water balance approach, may also be applied for determining the aquifer recharge when other quantities of water balance are well-known with high accuracy [16]. However, the outputs of the model are limited in the absence of model calibration and model validation due to the considerable uncertainty in input parameters. In view of this, for simulating spatial and temporal long-term groundwater recharge, empirical models have proved their authenticity in macroscale catchments [11,17]. Geohydrological models have also been used in regional estimation of groundwater recharge [18], where their applications include zoning via meteorological landscape and soil attributes. Recharge has been considered to be deeply spatially and temporally variable and not a linear function of climate, soil group types, land-use types, and other parameters of a drainage area [19]. Against this background, techniques have been developed in order to combine the spatial and temporal patterns of recharge [20,21]. Therefore, the climate and land-use effects on water budget quantities [22,23] can be spatially and temporally assessed [17] at the catchment or regional scale.

For the determination of groundwater recharge based on an empirical water balance model, a large variety of climate and catchment data (land cover, soils groups, geologic data, and topography) is required. The variety of the input data, e.g., land use, soils, and topography with various spatial and temporal ranges, together with various hydrological data requirements, complicates the interpretation of the outputs in terms of dominating factors. To generate common tendencies, the complexity and the dimensional space in data may be reduced in order to make interpretation easier [24]. Against this background, principal component analysis (PCA) can be applied, defined as a multivariate statistical method which applies orthogonal alteration to turn a set of the most correlated parameters into a set of linearly noncorrelated parameters, known as principal components or eigenvectors, while keeping the maximum amount of total variance in the dataset [24]. Every component is taken from a set of interrelated variables, which are affected by identical process, e.g., recharge and catchment data.

The main objective of this study is to quantify the groundwater recharge and to assess the recharge-controlling factors. In this study, a newly developed grid-based water balance model, mGROWA [11], was applied, due to its accuracy in simulating the spatially distributed water budget components, including groundwater recharge at a large scale, to determine the actual evapotranspiration, total runoff, direct runoff, and recharge in the Ergene catchment. The results of the model were interpreted by developing a methodology for assessing the most dominant factors influencing the groundwater recharge. For this purpose, principal component analysis (PCA) was applied to filter the data; accordingly, important independent basin factors presenting high variance in groundwater recharge were detected. The study was split into five sequential tasks: (i) to create a regular and reliable GIS input database made up of climate data (e.g., precipitation), land cover data, type of the soil, and topographical and geological data at the catchment scale; (ii) to perform, at a macroscale catchment, water budget analysis with the purpose of determining the renewable water resources, e.g., total runoff; (iii) to quantify predominating runoff components (direct runoff and groundwater recharge) by applying the base flow index (BFI) technique; (iv) to calibrate and validate the mGROWA model outputs by referring to a study area database of recorded runoff from different gauging stations; and (v) finally, to apply statistical methods to determine the dominant catchment and climate factors influencing groundwater recharge. To our knowledge, this is the first research to determine the temporally and spatially distributed groundwater recharge for sustainable groundwater reservoir management in the Ergene catchment.

2. Materials and Methods

2.1. Study Area

The Ergene catchment is situated in the middle of the Thrace zone, at the north-west part of Turkey. The catchment is bordered by the north Marmara catchment, Meriç catchment, and Bulgaria, and is located approximately between 40° 48’ 0” and 42° 09’ 0” N (latitude) and 26° 19’ 0” and 28° 11’ 0” E (longitude) (Figure 1). It is spread over 1.4% of Turkey and covers an area of about 11,021.5 km². Its topography reveals that the minimum and maximum elevations are 2.3 m and 1020.6 m above sea level, respectively. The main river of the catchment flows from the West part in the Istranca Mountains to the East part of the basin with 265 km of the length and the mean annual flow rate of 28.7 m³/s. The Ergene River and its affluent streams are significant sources of water in the north-western part of Turkey (Figure 1). The Ergene river catchment is dominated by terrestrial climate, summers are hot and arid in the north, whereas the winters are cold. In the south part, it is characterized by a Mediterranean climate where the summers are hot and dry and winters are warm and rainy. During the year, January and February are the coldest months, with an average minimum temperature of 4 degrees Celsius, whereas June and July are the hottest months, with an average maximum temperature of 20 degrees Celsius. The average annual temperature of the Ergene river basin is around 13 degrees Celsius. The mean annual long-term precipitation is about 590 mm, the highest precipitation amounts occur in November and the minimum in August. In the south, annual precipitation distribution increases from 650 mm/year to 730 mm/year, whereas in the northern part, the precipitation ranges from 600 mm/year to 900 mm/year [25]. In the Ergene river basin, since 1991, the precipitation information has been collected and no remarkable trend was found in precipitation [26,27]. The annual average evapotranspiration ranges from 450 to 950 mm. Hydrogeological units in the Ergene river basin area mainly composed by metamorphites (gneiss, schist, and marble) and clastic rocks [5,25]. The southeastern part is dominated by basalt fields, whereas alluvium exists in the valley floors, and tectonic lines cover the northern and northeast parts of the investigated area [5,25].

Figure 1 shows the Ergene river catchment geographic context, subcatchments, the location of available river gauging stations and monitoring wells. In Figure 1, the investigated sub-catchments, including Babaeski (488.08 km²), Inanlı (0.99 km²), Lulerburgaz (2697.16 km²), soğucak (71.45 km²), Çayırderekoy (54.48 km²), Ayvacık (1.3 km²), Payralı (90.38 km²), K.yoncak (3.37 km²), Inecık (100 km²), and Üzunköprü (10,026 km²), are presented. These subcatchments were delineated according to the available stream gauging stations with the ArcGIS watershed tool. The 32 monitoring wells existing in the investigated area are also presented in Figure 1. In this study, the groundwater depth data were collected from these wells for the period 1991–2010, where at least one record per month was available for each well.

2.2. Water Balance Modelling Based on the mGROWA Model

mGROWA is a grid-based hydrological water balance model developed for high resolution, long-term simulation of direct runoff, evapotranspiration, and groundwater recharge, and is applicable in macroscale catchments and on a regional scale [11]. mGROWA, in the first run, physically determines the total runoff and actual evapotranspiration in daily steps and the dynamics of soil water in the root zone (Equation (1)).

$$\Delta S=p+q_{cr}-e{t}_a-q_t$$

(1)

where S is the storage or soil moisture content in the root zone (mm) of areas covered with vegetation or the water stored on impervious surfaces in the case of urban areas. p stands for the precipitation in a grid cell (mm/day), q_cr is the capillary rise from the groundwater, et_a is the actual evapotranspiration (mm/day) determined by using Equation (2), and q_t is the total runoff (mm/day).

$$e{t}_a=e{t}_0\cdot K_{LN}\cdot f\left(\beta ,\gamma \right)\cdot f\left(s\right)$$

(2)

where et₀ is the grass reference evapotranspiration (mm/T), which, in this study, was estimated using a modified Hargreaves and Samani approach [28], K_LN stands for crop coefficiency and corrects et₀ respective of land-use. From this perspective, it considers various species of vegetation with properties different from the reference grass area. f(β,γ) is a terrain characteristic used to modify et_a according to the slope γ and exposition β, and f(s) stands for a storage function.

$$f\left(\beta ,\gamma \right)=\gamma \left[1.605\cdot {10}^{-2}\cdot \sin \left(\beta -90\right)-2.5\cdot {10}^{-4}\right]+1$$

(3)

f(β,γ) was obtained by applying Equation (3) [29].

Figure 2 presents the aspect and slope of the Ergene river catchment. It can be seen that the minimum slope in the study area is 0%, whereas the maximum slope is 66.5% in the northern part of the catchment. The arithmetic mean slope was observed to be 5% with a standard deviation of 4.3%. The aspect reflects the structure of the land surface in the investigated area. The spatial model’s resolution is arbitrary and can be freely chosen. For the Ergene catchment, a 100-m grid is considered suitable, resulting in a total number of 1,036,861 grid cells. The water balance quantities were calculated individually for each grid cell and only the vertical motion is considered as the grid cells are not interconnected; therefore, the lateral flow is not counted in mGROWA.

In the second step, mGROWA empirically divides the total runoff into direct runoff and groundwater recharge based on monthly time steps as the groundwater management does not require a higher temporal resolution [11]. The separation of total runoff into dominant components (direct runoff and groundwater recharge) is based on the theory of base flow indices (BFI) as described in Kunkel and Wendland, (2002) (Equation (4)) [16].

$$q_t=BFI\cdot q_t+\left(1-BFI\right)\cdot q_t=q_r+q_d$$

(4)

where q_d represents direct runoff and q_r is the groundwater recharge.

The attribution of area-differentiated baseflow indices by mGROWA is based on BFI from literature depending on the type of rocks (unconsolidated and hard rock types).

Table 1. Base flow index (BFI) values of the unconsolidated rock areas according to Dörhöfer and Josopait (1980) [31], Hennings (2000) [32], and Wessolek and Facklam (1997) [33].

Degree of Sealing	Groundwater Depth	Waterlogging Tendency	Slope	Baseflow Indices
	<2 m	No water logging	<1%	1
	1.3–2 m	1 (very low)	1–3.5%	0.9
I (10–45%)				0.82
			3.5–7%	0.67
			7–10%	0.59
	0.8–1.3 m	2 (low)	10–13%	0.5
	0.4–0.8 m	3–4 (medium – high)	13–15%	0.44
	<0.4 m	5 (very high)	>15%	0.4
II (45–75%)				0.33
III (75–90%)				0.28
IV (>90%)				0.2

shows the baseflow indices values used in mGROWA for unconsolidated rock areas. As can be seen in Table 1, the most important attributes for estimating baseflow indices are groundwater table, water logging tendency of the soil, the degree of sealing, and the area-related slope gradient. In the solid rock areas, the baseflow index values presented in Table 1, may not be applied due to the influence of geological subsoil conditions that predominate the soil properties in quantifying the runoff fractions [30].

Table 2. The classification of the hydrogeological properties of the hard rocks in North Rhine-Westphalia and the associated BFI values obtained by calibration [20].

Hydrogeological Class	Permeability	Hydraulic Conductivity	Baseflow Indices
I	Very high	>10−2 m/sec	0.9
II	High	>10−3–10−2 m/sec	0.6
III	Medium	>10−4–10−3 m/sec	0.57
IV	Moderate	>10−5–10−4 m/sec	0.3
V	Low	>10−7–10−5 m/sec	0.29
VI	Very low	>10−9–10−7 m/sec	0.18
VII	Extremely low	<10−9 m/sec	0.12

presents the categories of geological characteristics of hard rocks, such as the permeability, as stated by Bogena et al. (2003) [34], and aquifer typologies initiated by Wendland et al. (2007) [35]. The baseflow indices approach is established by assuming that aquifers possessing similar petrographic characteristics exhibit similar or comparable hydrodynamic and hydrologic conditions regardless of their regional occurrence worldwide. Then, baseflow indices are given for single model grid cells of a groundwater recharge model depending the geological information. This practice has had successful application in Germany [11], Turkey [36], Slovenia [37], Greece [38], and Cameroon [39], yielding realistic groundwater recharge values.

Figure 3 shows the base flow index (BFI) methodology applied in mGROWA for simulating the groundwater recharge in the Ergene catchment. In the determination of groundwater recharge, the hierarchical method is split into five tasks. Firstly, in every grid cell, it is supposed that recharge is negligible for an impervious surface (1). Subsequently, it checks the remaining grid cells for the existence of significant artificial drainage, such as channel or ditch drainage, and its corresponding base flow index is determined (2). On the condition that there is no impervious surface or artificial drainage in the investigated area, the analysis of the type of the rocks (consolidated or unconsolidated rock) is done (3). For consolidated rock, attributes related to the hydrogeology of the case study are applied in BFI estimation (4). For unconsolidated rock regions, the BFI is estimated by taking into account the groundwater table levels, the slope, and water logging tendency (5). Accordingly, the impervious surface, e.g., urban areas, display a total runoff which is equal to the direct runoff. Under that surface, the recharge is equivalent to zero as on the impervious layer; no penetration of the water occurs. In the pervious surface, e.g., vegetation cover, the groundwater table level plays an important role in separation. In a shallow zone, where the influence of capillary rise is observed, the net groundwater recharge will be determined. Where the influence of capillary rise is not considered, the rock types control the recharge. Depending on the types of rocks, the BFI can be calculated [16].

2.3. Baseflow Determination

There are only ten streamflow observation stations available with a long hydrograph record for the period of 1970–2014. The related gauging stations and subcatchments were delineated by means of the hydrology tool in Arc-GIS on the basis of the digital elevation model and river grid map of the Ergene catchment (Figure 1). For the validation assessment of mGROWA outputs, the daily mean of total runoff at selected gauging stations was evaluated. The ten streamflow stations constitute different areas within the basin, from 1 to 10,026 km², and for a great variety of climatological, geological, and site soil conditions. The observed mean long-term total runoffs were applied to check the validity of the simulated mean long-term total runoff. Groundwater recharge can only indirectly be determined from the recorded streamflow by applying relevant separation methods. The base-flow for selected subcatchments was separated from streamflow based on monthly low water flow rates (MoLR) [40]. Wundt (1958) [40] demonstrated that in unconsolidated aquifers, for prolonged series of years (20 years), the mean MoLR values are good estimates of groundwater recharge. This method is referred to as the MoMLR method, where MoMLR stands for Monthly Mean Low water discharge for prolonged series of years. In hard rock aquifer areas, however, the MoLR-values do not fit to the average groundwater recharge, as MoLR includes considerable components of direct runoff; e.g., direct runoff and interflow, considered also by the MoMLR method. Therefore, Kille (1970) [41] and Demuth (1993) [42] developed a suitable hydrograph separation method applicable to consolidated aquifers, to permit a reduction of MoMLR values by the interflow quantity.

After determining the regional averages of the base flow index values for the investigated basin areas based on river discharge, it was followed by the disaggregation and the transfer base flow indexes to areas where there is no available suitable gauge data. Thus, there was a consideration of parameters and runoff-effective regional features where the hierarchical approach in Figure 3 was taken into account, in which the value of a site condition is taken into account to be exclusively estimating the base flow index value. In the procedure of iterative calibration [36], BFI was calculated by adopting Kille’s (1970) techniques [41], and various characteristics (degree of sealing, groundwater depth, water logging tendency, slope, and permeability) were identified. Among the total gauging stations with possessing data, half of the drainage points were used for calibration and the other half were used for validation of groundwater recharge [43]. Accordingly, a set of base flow index values is established and guides the most favorable fit for the investigated area.

2.4. Statistical Analysis

The comparisons of simulated and recorded water balance quantities in the Ergene catchment were assessed using the Nash–Sutcliffe index [44] to spatially evaluate divergences in simulated groundwater recharge and total runoff with recorded hydrographs in subcatchments possessing sufficient stream flow observations. The Nash–Sutcliffe index, also referred to as the efficiency index (E_f), has been widely applied and is a valid statistic indicator for evaluating the goodness of fit of a model [45]. E_f may vary in the range of −∞ and 1. The acceptable range of agreement varies from 0 to 1, a negative value reveals disagreement of the performance of the hydrological models, and a value of 1 shows excellent performance. Herrmann et al. (2015) [11] presented a different concept for simultaneously applying E_f to the all considered stream flow stations and their corresponding drainage areas within a catchment by applying Equation (5).

$$E_f=1-\frac{{{\displaystyle \sum }}_{i=1}^n(A_i\cdot {\left(Q_{observed,i}-Q_{simulated, i}\right)}^2)}{{{\displaystyle \sum }}_{i=1}^n(A_i\cdot {\left(Q_{observed,i}-Q_{observed,A_V}\right)}^2)}$$

(5)

where Q_observed stands for the recorded stream flow value per unit of area, Q_simulated the related modeled stream flow, $Q_{observed,A_V}$ the recorded overall flow rate per unit of all sub-basins taken into account. A is the specific basin area and i the basin indicator.

In this study, a multivariate technique called principal component analysis (PCA) was adopted through the XL-STAT software. PCA alters an original set of dependent variables into an equivalent number of noncorrelated factors, referred to as principal components (PCs). In this study, PCA was adopted to reduce the amount of redundant data and to determine the dominant variables that show most of the variance of groundwater recharge, thus the factors most importantly influencing the spatial variation of recharge. Its purpose is to withdraw the important information from the dataset. This approach squeezes the data by maintaining only the most significant information, facilitating the explanation, and investigating the form and constituents of the data [46]. A varimax rotation method was conducted in order to increase the variation between the parameters for every component [47].

In the first step, the strength of existing correlations between groundwater recharge and catchment attributes were determined by the Pearson product–moment correlation coefficient [48]. For the matrix of interrelated variables, the linear relationship among the catchment attribute was tested by measuring the correlation between parameters in the range of [−1,1]. Greater coefficients show the importance of the correlation and high interdependence among the variables. A positive coefficient represents similarity and harmony among the interrelated variables, and a negative sign reflects the opposite [49].

In the second step, the variables that displayed significant correlation with the recharge were analyzed by applying PCA. The first PC reflects a great portion of the variance in the original information. The PCs possessing an eigenvalue greater than one and containing a high factor loading (>0.7) were considered as the variables that control the spatial variation of recharge [12,17]. Following PCA, the highly weighted variable in a principal component was determined based on the level of correlation coefficients and their sums [48]. This approach minimized the number of variables to produce a small data set that highly influenced recharge.

2.5. Input Data

All the climatic data were recorded from 41 stations located in the Ergene catchment and were provided by the Turkish General Directorate of State Hydraulic Works (DSI). The database consists of daily precipitation and daily extreme temperature data. The latter data was used to determine the reference grass evapotranspiration level at corresponding stations by applying Hargreaves’ method [50]. This method was established by only referring to the mean air temperature of a region, temperature extremes, and extraterrestrial radiation. It is approved by the Food and Agriculture Organization of the United Nations (FAO) for its reliability and accuracy in estimating reference evapotranspiration where there are gaps in climatic data; e.g., in relative humidity, net solar radiation, and wind speed [50]. However, in humid regions, the Hargreaves equation requires calibration [51]. Cobaner et al. (2017) [28] proposed a modified Hargreaves equation to improve its estimation throughout Turkey. The dataset was produced as digital datasets for the hydrologic period of 1991 to 2010 in a gridded pattern for the study area.

The data was used to produce catchment-wide 100 × 100 m data by applying geographic information system (GIS) interpolation techniques saved in an American standard code for information interchange (ASCII) text grid file applying the unique value. Precipitation and grass reference evapotranspiration values were interpolated using inverse distance weighted techniques (IDW) [52] and their values were allocated to every single grid (100 × 100 m) of the catchment in ArcGIS. This method was applied to produce the area distribution of the precipitation and grass reference evapotranspiration in the Ergene catchment. The climate input data was developed spatially and temporally, such that long-term mean precipitation and grass reference evapotranspiration for the period of 1991–2010 over the study area was approximated to 591.29 ± 66.9 mm/ year and 1035.26 ± 34.46 mm/year, respectively. The more standard deviations signify that the spatial variability of climatic data is significant over the catchment. Then, the spatially distributed climatic data were used in mGROWA for the entire study catchment. The spatially distributed daily mean climatic data was aggregated into monthly values and the grids of monthly sums of precipitation and grass reference evapotranspiration in mm were prepared as model input.

Table 3. Database of the mGROWA water balance model for the Ergene catchment.

Data	Data Base	Scale/Spatial Resolutions	Data Source
Climate data	Precipitation (1991 to 2010) Temperature (1991 to 2010)	100 × 100 m	General Directorate of State Hydraulic Works (DSI)
Soil data	Main soil groups	1/25,000	National soil database by General Directorate of State Hydraulic Works (DSI)
	The depth of the groundwater table, Perching water influence Root depth, effective field capacity	Derived based on pedo-transfer functions	National soil database by General Directorate of State Hydraulic Works (DSI) and fields studies
	Capillary rise from the groundwater table
Land cover	Land use types Percentage imperviousness	1/25,000	National soil database by General Directorate of State Hydraulic Works (DSI)
Hydrogeology	Geological map	1/500,000	General Directorate of State Hydraulic Works (DSI) and fields studies
Relief	Digital elevation model SRTM Ground surface slope Ground surface exposition	30 × 30 m	SRTM/X-SAR by National imaging and mapping agency
Runoff		Station data Monthly and daily resolution	General Directorate of State Hydraulic Works (DSI) and fields studies

shows the prepared dataset of mGROWA for the study area.

The land use information was required in order to estimate the real vegetation evapotranspiration and was collected from the digital landscape model. The concept of CORINE (coordination of information on the environment) land cover contains 44 land cover classes, of which 26 occur in the study area (Figure 4).

Figure 4 presents the land use types of the Ergene river catchment for the year 2012. In Figure 4, one can see that the majority of the Ergene catchment is covered by vegetation, i.e., nonirrigated arable land of code 211 (approx. 54.9% of the total study area), followed by permanently irrigated land of code 212 (approx. 8%), land principally occupied by agriculture of code 243 (approx. 7.1%), transitional woodland-shrub of code 324 (6.8%), broard-leaved forest of code 311 (6%), pastures of code 231 (3.9%), natural grasslands of code 321 (3.2%), rice fields of code 213 (2.6%), and discontinuous urban land (112) makes up 1.7%. The share of other categories exhibits low percentages (below 1%). For Ergene water balance modeling, the information of imperviousness, rooting depth of plants, and land-use-specific actual evapotranspiration factors were allocated to the Ergene digital landscape model correspondingly.

Figure 5 shows the soil groups in the Ergene river catchment. As can been observed from Figure 5, in the Ergene catchment, eight soil groups exist. In the north of the study area, brown forest soil without lime mostly dominates. The central part of the Ergene catchment is covered by brown soils without lime and vertisol groups. The river beds and its surrounds are characterized by the existence of alluvial soils. The southern part of the study area is covered by brown forest soils with and/or without lime. The urban soils are distributed in small quantities over the catchment. The soil groups mainly occurring in the study area are brown soils without lime, brown forest soils without lime, and vertisols, respectively contributing 33.8, 25.2, and 23.6%, followed by alluvial soils (9.2%), brown forest soils (5.8%), urban soils (2.1%), and rendzinas soil (0.02%). The parameterization of Ergene soil profiles has been done by applying pedo-transfer functions. The rooting depth and effective field capacity were derived to present the storing and releasing capacity of the Ergene soils.

3. Results

3.1. Simulation of Water Balance Quantities

The long-term annual mean of spatially distributed water budget components calculated with the mGROWA model over the study area for the period 1991–2010 are presented in this section.

Figure 6 presents the water budget components: the actual evapotranspiration (ET_a), runoffs, and recharge. The runoff maps reveal the individual runoff quantities produced in 100 × 100 m grid cells; therefore, the Ergene River flow rates passing through the investigated area are not included. The actual evapotranspiration is predominantly affected by vegetation cover, soil, and climate. It can be noticed that the area-distributed ET_a presents considerable variability in land cover across the study area. The lowest ET_a (less than 300 mm/year) was observed where the impervious areas dominate, especially in urbanized areas. A high ET_a in the range of 400–450 mm/year appears in the highest altitudes of the northern part where the land is covered by forest. In the low altitude land, ET_a is also influenced by the land-use types. In the nonirrigated arable land (dominating in the catchment (55%)), the ET_a varies in the range of 350–400 mm/year. The arithmetic mean value of the real evapotranspiration is 450 mm/year over the study area.

Regarding the runoffs (Figure 6), because of the high level of imperviousness, the urbanized areas display the highest total runoff values, with direct runoff levels in the ranges of 350–400 mm/year and 200–300 mm/year, respectively. The low permeability of urbanized areas promotes runoff. The total runoff levels steadily decrease in the high-altitude land located in the northern part due to the forest cover which generates high actual evapotranspiration (less than 150 mm/year). Accordingly, the highland areas in the north-eastern part display total runoff values less than 100 mm/year, whereas in the north-western part of study area, values in the interval of 100–150 mm/year prevail. The generated direct runoff shows similar spatial patterns to the total runoff and both are highly influenced by the land-use type. Across the Ergene catchment (approx. 1,091,645 grid cells), the arithmetic averages for total runoff and direct runoff are 200 mm/year and 85 mm/year, respectively.

Regarding simulated long-term groundwater recharge in Figure 6, the Ergene catchment can be mainly characterized into three zones: high, moderate, and low groundwater recharge rates depending on the zone variety parameters, e.g., soil, geology, the ability of the aquifer to store water, hydraulic conductivities. The influence of hydrogeology is obvious, i.e., the alluvial plains in the lowland of the catchment prevail a high groundwater recharge levels greater than 350 mm/year, because of high infiltration rates in such unconsolidated sand aquifers. Besides, the southern part receives high precipitation, thus, some grids located in alluvial soils display the maximum groundwater recharge values. In the moderate zone, groundwater recharge displays mainly values in the interval of 100–200 mm/year. In the northern zone of the catchment, the groundwater recharge rates are generally diminishing, down to the lowest values (less than 50 mm/year). The calculated long-term average groundwater recharge level over the study area is 95 mm/year.

The results for the average monthly groundwater recharge for the period of 1991–2010 are presented as a map in Figure 7. It can be observed that the groundwater recharge rates start to occur in November and end in April. The highest rates were observed in January. In April, the groundwater recharge starts to decrease as the vegetation time begins.

In Figure 8 and Figure 9, the average annual water budget quantities—actual evapotranspiration, direct runoff, and groundwater recharge—are displayed in terms of existing CORINE land cover classes and soils group types in the Ergene catchment. Their corresponding coverage areas in the catchment are also shown.

In Figure 8, the groundwater recharge shows great variation depending on the land-use. The artificial surface, represented by CORINE_2012 codes 111, 112, 121, and 122 (continuous-discontinuous urban fabric and industrial areas) carry a low groundwater recharge value due to the high level of imperviousness characterizing these types of land-use. Consequently, these land-use types generate high runoff. Descriptions of the codes are presented in Figure 4. Artificial land covers 3% of the area of the Ergene catchment. The higher levels of groundwater recharge can be observed in land use types such as mineral extraction sites, dump sites, and constructions sites (represented by Corine codes 131, 132, and 133, respectively) which only comprise 0.3% of the study area. These land-use types (131, 132, and 133) affect the natural hydrogeology quite profoundly, where the permeability may increase up to 100 times [53,54]. For agricultural areas (211, 212, 213, 221, 222, 231, 242, and 243), the actual evapotranspiration dominates other water balance quantities, which is also the case in the forest and seminatural cover areas (311, 312, 313, 321, 323, 324, and 334). These types of land-use cover displayed a moderate average level of recharge. The exception was observed in bunt cover areas (334), where groundwater recharge (240 mm/year) is higher than actual evapotranspiration.

In Figure 9, the groundwater recharge is shown to be dependent on soil group types in the study area, which is dominated by alluvial soils. This type of soil presents a high ability of water infiltration to replenish the groundwater aquifer. A low average groundwater recharge of 50.3 mm/year is observed for the urban soil type. Urban soil is a type of anthropic soils, widely formed by human activities, in urban areas.

Figure 10 illustrates the descriptive statistics of base flow obtained from the Wundt and Kille separation techniques. The highest median was noticed in Soğucak, Çayirderekoy, Ayvacık, and Poyralı subcatchments of 90.4, 77, 73.5, and 49 mm/year, respectively. The mentioned subcatchments are found in the northern part characterized by high topographic elevation and slope (Figure 2). The relief significantly affects the baseflow and recharge. The increase in high relief, e.g., slope angle, amplifies the level of groundwater transmission [55]. Other subcatchments, especially Luleburgaz, Babaeski, and Inanlı, displayed the lowest minimum median in estimated base flow with 33, 37, and 38 mm/year respectively. These subcatchments are located in areas of medium slope and medium elevation. On the other hand, the low value of base flow in these subcatchments may be justified by their location on land covered by agricultural land (73%) and are frequently watered from groundwater. A portion of urban areas existing in these subcatchments may also negatively influence base flow due to the change in the natural system’s surface and subsurface capacity of to drain water. As a result, the interaction between groundwater recharge and rivers is deteriorated [56]. In the Ergene catchment, it can be concluded that catchment characteristics (e.g., topography, slope land-use types) and the interaction between surface and subsurface water manage the recharge-base flow interconnection.

3.2. Model Evaluation

mGROWA cannot be characterized as a runoff concentration model; therefore, its validation does not rely on comparing daily simulated and observed water balance quantities [11]. However, validation of the mGROWA model outputs is done by comparing the long-term observed and simulated runoffs and recharges. Accordingly, Figure 11 shows the scatter plots with a 1:1 line, which presents the comparisons between modeled total runoff and groundwater recharge with their corresponding observations for the hydrological period 1991–2010. It can be observed from the left of Figure 11 that the comparison displays a quantitative Nash–Sutcliffe statistic (E_f) of 0.94 for observed-simulated total runoff, which indicates a very good performance of the hydrological model. On the right side of Figure 11, the good performance of the mGROWA model (E_f = 0.76) for calculated groundwater recharge and base flow obtained from the Wundt–Kille separation techniques was achieved and taken as observed data. As stated by Moriasi et al. (2007) [57], a Nash–Sutcliffe efficiency value greater than 0.75 and 0.65 shows a very good and a good performance of the hydrological model in an individual basin, respectively. As can be observed, the deviations between observed and simulated total runoff and groundwater recharge for nearly all considered subcatchments are below ±20%. For only three of the sub-catchments considered are the deviations between the recorded and modeled groundwater recharge levels higher than 20%. The inaccuracy in simulated values and observations may be caused by scale discrepancy of the empirical model [20]. In this comparison between the observed and the model outputs, there are obvious differences. These inaccuracies can be caused by unreliability in the complexity of the model-data environment, such as the modelling approach, errors occurring in observation, uncertainties occurring in soils groups, land-use types, topography, and climatic data. In addition, the estimated groundwater divides in some considered subcatchments might not be representative for real corresponding river drainage area due to artificial channels, irrigation practices, dams, and pumped groundwater wells that can perturb the natural water budget of a catchment. Besides, the fact of assuming that recharge is equivalent to base flow implies the equivalence between the groundwater shed and the surface watershed, where the groundwater cannot discharge across boundaries. Although groundwater can run into an aquifer, most of the time it overlaps many basins, and within basins, groundwater and surface water boundaries frequently differ from each other, and this may change in time. Therefore, identifying the groundwater shed area that contributes to a river is a hard work. Usually, a dissimilarity between the observed base flow in the river and the estimated recharge may be caused by a lack of consideration of the lateral fluxes in such models that consider only vertical fluxes, and the effect of scale.

3.3. Determination of Dominant Factors

Prior to principal component analysis, the Pearson correlation test was applied to indicate the degree of linear association between catchment characteristics (soils, land use, area, topography, and climate) and the groundwater recharge. Twenty catchment attributes were assembled for processing by means of Arc-GIS tools.

Table 4. Catchment attributes, units, and their Pearson correlation with groundwater recharge.

Catchment Attributes	Abbreviation	Units	Correlation with the Groundwater Recharge
Soil types
Alluvial soil	AS	%	0.671
Brown forest soils	BFS	%	−0.293
Brown forest soils without lime	BFSWL	%	−0.752
Rendzinas	R	%	−0.277
Brown soils without lime	BSWL	%	−0.020
Vertisol	V	%	0.302
Urban soils	US	%	−0.062
Land cover characteristics
Artificial surface Land	ASL	%	0.130
Agriculture area	AA	%	0.676
Forest	F	%	−0.659
Semi natural area	ANA	%	−0.300
Wetland and water bodies	WWB	%	−0.072
Area and topography
Basin drainage area	BDA	Km2	0.541
Basin slope	BS	%	−0.058
Elevation	E	M	−0.469
Depth to water table	DWT	M	0.185
Climate
Precipitation	P	mm/year	0.817
Actual evapotranspiration	AE	mm/year	−0.686

Note: values in bold represent correlation with a significance level alpha = 0.05.

shows the catchment characteristics, symbols, units, and their correlation with groundwater recharge. The catchment characteristics are mostly presented as a percentage of their corresponding areas relative to the total area of the Ergene catchment. In Table 4, the correlation between the catchments’ characteristics and groundwater recharge are presented. The correlation assessment in

Table 5. Correlation analysis of catchment characteristics.

Variables	GR	BSWL	AS	BFS	V	BFSWL	US	R	AS	AA	F	ANA	WWB	BDA	BS	P	AE	DWT	E
GR	1.00
BSWL	−0.02	1.00
AS	0.67	0.21	1.00
BFS	−0.29	−0.61	0.11	1.00
V	0.30	−0.04	−0.16	−0.43	1.00
BFSWL	−0.75	−0.23	−0.30	−0.12	−0.53	1.00
US	−0.06	0.20	0.47	−0.21	−0.25	0.23	1.00
R	−0.28	−0.30	−0.42	0.41	−0.20	0.07	−0.16	1.00
AS	0.13	0.43	0.55	−0.47	0.15	−0.13	0.62	−0.23	1.00
AA	0.68	0.41	0.33	−0.34	0.73	−0.85	−0.08	−0.46	0.34	1.00
F	−0.66	−0.44	−0.37	0.38	−0.72	0.83	0.02	0.47	−0.42	−1.00	1.00
ANA	−0.30	0.20	0.18	−0.33	−0.16	0.36	−0.03	−0.16	0.18	−0.24	0.20	1.00
WWB	−0.07	0.29	−0.15	−0.57	0.45	−0.08	−0.27	−0.32	−0.16	0.30	−0.29	0.57	1.00
BDA	0.54	0.24	0.53	−0.28	0.14	−0.15	0.28	−0.15	0.37	0.34	−0.36	−0.03	0.06	1.00
BS	−0.06	−0.39	0.07	0.49	−0.56	0.43	−0.25	−0.17	−0.59	−0.54	0.58	0.35	0.08	−0.24	1.00
P	0.82	−0.22	0.51	0.59	−0.01	−0.53	−0.35	−0.25	−0.33	0.34	−0.29	−0.17	−0.07	0.19	0.44	1.00
AE	−0.69	−0.42	−0.60	0.25	−0.43	0.69	−0.30	0.30	−0.72	−0.80	0.83	0.25	0.11	−0.57	0.64	−0.19	1.00
DWT	0.19	0.31	0.12	−0.52	0.70	−0.50	0.06	−0.13	0.68	0.64	−0.69	0.15	0.28	0.12	−0.69	−0.26	−0.68	1.00
E	−0.47	−0.58	−0.39	0.45	−0.57	0.72	0.02	0.40	−0.45	−0.85	0.87	−0.24	−0.55	−0.23	0.43	−0.17	0.68	−0.75	1.00

Note: values in bold represent correlation with a significance level alpha = 0.05.

indicates that among the variables, significant correlations were observed. Consequently, some variables are absolutely redundant. In Table 4, only six catchment characteristics were found to be significantly correlated with groundwater recharge. The climate parameters appeared to influence the area distribution of recharge. The highest positive correlation was obtained by precipitation (0.817), and this was found to be a natural principle in many hydrogeological studies [12].

The actual evapotranspiration showed a negative correlation with recharge. This means that high actual evapotranspiration will deny the aquifer maximum groundwater recharge and low actual evapotranspiration will cause an increase in groundwater recharge. The same finding has been previously presented [58]. It was followed by Brown forest soils without lime, which was detected to be inversely correlated with groundwater recharge. Two types of soil appear to be strong in a relationship with recharge. Alluvial soil presents a strong positive correlation with recharge. The Ergene basin’s alluvial aquifers are formed principally by sandy soil which possesses intermediate to high hydraulic conductivity, thus has high infiltration capacity, leading to an increase of groundwater recharge. Brown forest soils without lime appear to correlate negatively with recharge. The most important constituent of this type of soil is the clay which has low permeability. Two of the five land-use types indicated a significant correlation with recharge. Agriculture land and forests were found to have significant interdependency with recharge. The agglomeration of the forest biomass is likely to increase the evapotranspiration and decrease the total runoff, thus negatively affecting the amount of water that replenishes the aquifer [59]. The area of the catchment and its topography did not show any significant correlation with the groundwater recharge. Principal component analysis (PCA) was applied to the variables most closely correlated with groundwater recharge. AS, BFSWL, AA, F, P, and AE displayed significantly interdependency with the recharge.

Table 6. Principal component analysis (PCA) of catchment attributes that are most correlated with groundwater recharge.

Statistical Procedure	Indicators	PC1	PC2
PC Analysis	Eigenvalue	3.946	1.090
	Variability (%)	65.774	18.174
	Cumulative (%)	65.774	83.948
Variables
AS		0.582	0.621
BFSWL		−0.894	0.087
AA		0.939	−0.277
F		−0.941	0.289
P		0.516	0.719
AE		−0.880	0.139

Note: bold type indicates variables with a loading factor greater than 0.70.

shows the principal component analysis of catchment factors displaying a high correlation with recharge. Eigenvalues greater than one were considered as significant in this study (Table 6). The spatial variance of the recharge for the two considered principal components was 84%. Percentage of forest land cover yielded the largest loading factor, with agricultural land cover coming second, followed by brown forest soils without lime and actual evapotranspiration, respectively, in PC1. The precipitation impacts the recharge more in PC2.

Catchment factors displaying high positive weighting on PC1 include alluvial soil (AS), agriculture area (AA), and precipitation (P) (Table 6). Precipitation showed a strong positive correlation with recharge since precipitation is the main source of groundwater recharge, as confirmed by many studies [17]. Alluvial soils increase the amount of recharge, as this soil type is associated with high rates of sand which allows infiltration into alluvial aquifers [60]. Therefore, it is plausible that alluvial soils are related to an increase in recharge. The positive correlation between recharge and agricultural land use may increase the recharge depending on its patterns of infiltration and evapotranspiration; nevertheless, without knowing the details of these land covers in the Ergene catchment, the direct mechanism cannot be ascertained. Although these interrelations identify the catchment factors which can control the spatial patterns in recharge, further study is recommended to determine the exact processes by which these factors contribute to the increase of recharge across the catchment.

The PCA outputs show that recharge is lower in catchments where brown forest soils without lime, forest, and actual evapotranspiration dominate in Table 6. Brown forest soils without lime were found to have a considerable rate of clay soils, which are highly related with evapotranspiration and runoff, which leads to decreased groundwater recharge. Actual evapotranspiration negatively affects the groundwater recharge by diminishing the available quantity of water for recharge. Forest cover is known to assist in infiltration and possesses a relatively high amount of evapotranspiration, thus low recharge. In accordance with PCA, the recharge variation is controlled, in order of significance, by vegetation land-use, soil group types, and climate. The great influence of vegetation on groundwater recharge was previously confirmed in a global assessment of groundwater recharge [61]. Globally, vegetated land provides the greatest amount of water input, such as precipitation and irrigation water, that becomes recharge [61]. The recharge in coarse soil is much greater than in fine textured soil, thus the movement of the water in the soil is highly controlled by soil texture and structure. This result is in agreement with previous studies [12]. The great influence of climate on recharge was also confirmed by Jan et al. (2007), in their study [62] about rainfall intensity and distribution impact on groundwater level changes, in which they found that the increase or decrease in level of groundwater is mainly a result of the presence or absence of precipitation, respectively. Thus, groundwater directly relies on accumulated precipitation.

4. Conclusions

The Ergene River catchment is one of the important basins in Turkey that are exposed to a high level of groundwater water exploitation due to industrial and agricultural activities. Thus, quantifying the groundwater recharge is important for its management. The spatial groundwater recharge is improperly considered in the modeling of groundwater recharge. However, hydrograph separation methods show spatial recharge distribution. With the development of GIS and the accessibility of spatial inputs, the recharge can be determined as the outcome of water budget modeling. The mGROWA model can, therefore, accurately determine the spatially distributed long-term groundwater recharge of the Ergene River catchment. The study outcomes show that the groundwater recharge at the scale of the Ergene catchment has significant variation, while in small-scale catchments, the variation is less noticeable. Mineral extraction sites, dumps, and construction sites display the highest recharge rates, followed by agriculture land. Regarding soil structure and texture, alluvial soils produce the highest groundwater recharge levels. Urban soil displays the lowest recharge rate. Accuracy is challenging in recharge determination. In this respect, a base flow assessment of ten available river gauging stations was used as a validation approach where the average base flow fluxes for a long period can be considered as a proxy for mean groundwater recharge rates. The assessment of groundwater recharge together with the catchment attributes, such as land-use, soil types, and climate variables, displayed negative and positive correlations. In accordance with PCA, the spatial groundwater recharge variation is controlled in order of significance by vegetation land-use, soil group types, and climate. To aid water resource management in Turkey, the approach developed in this study is suggested to be applied in different catchments across the country. The results may contribute to the establishment of a quantitative status of groundwater water reservoir exploitation and enhancement of the water resource management of the region.