Variability of Intrinsic Groundwater Vulnerability to Pollution in River Valley due to Groundwater Depth and Recharge Changes

Groundwater management can be effectively implemented by assessment of groundwater vulnerability to pollution. The research was carried out in the Vistula River valley (Poland) in an area of shallow groundwater occurrence, defined as a Groundwater-Dependent Ecosystem area. The goal of this study was to identify the average, maximum, and minimum depths of the groundwater table for variant assessment of groundwater vulnerability to contamination. The variants correspond to the average vulnerability, the vulnerability during hydrogeological drought, and the vulnerability during the flood risk period. Theoretical and effective weights of vulnerability parameters were calculated using the DRASTIC method. Vulnerability maps constructed for the various vulnerability variants and by using different parameter weights show the spatial variability of the individual vulnerability classes. Due to the specifics of this protected area, the expected dependency between vulnerability index and chloride concentrations in the monitoring points was not found. Comparison showed strong dependency of water chemistry and the value of recharge, lithology of aquifer, and unsaturated zone. The research results confirm the need for variant vulnerability assessment to protect against floods or predict the effects of climate change.


Introduction
Maps of groundwater vulnerability to contamination, which are the main tool in spatial planning and prediction of the effects of aquifer pollution, are compiled using many methods and procedures.Parameters and data taken for the vulnerability assessment are represented by average values and therefore are characterized by an "average" groundwater vulnerability to contamination.In the assessment of natural groundwater vulnerability to contamination, DRASTIC is the rating method commonly used.It is dedicated to areas of porous aquifers [1].In the DRASTIC method, like in other methods of groundwater vulnerability assessment, the depth to groundwater table (D) and the groundwater recharge by infiltration (R) depend on climate, weather, hydrological, and environmental changes and are the main parameters of vulnerability assessment, especially in shallow groundwater areas.The value of D is taken for the vulnerability assessment based on statistical data from long-term monitoring [2][3][4][5][6][7][8], and the value of R is based on averaged precipitation data, lithologies of the near-surface zone, or retention data [4,6,[9][10][11][12].The pollution vulnerability can be defined as 'the sensitivity of groundwater quality to an imposed contaminant load' depending only on the 'proper characteristic of the aquifer'.The vulnerability is relatively static, excluding some variation over time due, e.g., to the piezometric level changes [13][14][15][16].
The aim of the studies was the identification of characteristic groundwater levels in the Vistula River valley for the variant vulnerability assessment.During the studies, the variants of groundwater vulnerability to contamination were calculated by taking into account the average, maximum, and minimum values of D. These variants correspond to the average vulnerability, vulnerability in conditions of hydrogeological drought [17], and vulnerability in the conditions of floods in river valleys [18].The groundwater vulnerability assessment was carried out for three variants of infiltration: for the average groundwater levels from the period 1999-2013, and for the lowest and highest (average annual) values of groundwater levels.Sensitivity analysis and map removal parameter for all the variants, interpretation of vulnerability variations due to groundwater and recharge changes, and vulnerability validation with hydrochemical data usage were carried out.

Study Area
The analysis of vulnerability assessment changes was performed on a well-examined area located in the Middle Vistula Valley, Poland.It has been covered by legal protection since 1959 as the Kampinos National Park [19], Biosphere Reservoir [20] and Natura 2000-Kampinos Forest [21].The national park covers an area of 385 km 2 .The vulnerability assessment was carried out in the part of the Vistula Valley that covers 571 km 2 .Valuable plant communities whose existence is ensured by the shallow groundwater table cause areas to be classified as groundwater-dependent ecosystems (GDE), which cover approximately 62% of the park acreage [22].The main factor determining proper functioning of these GDEs is dependency: precipitation-infiltration-ecosystem type-surface water drainage.In the GDE areas, about 33% is covered by typical swamp vegetation, 6% is represented by wet communities, and 19% are dry communities, where the humidity of the aeration zone and the depth of the groundwater table are insufficient to maintain the hydrophilic vegetation [22].
In the Vistula Valley, Quaternary sediments are the collector of groundwater.The surface of tertiary loams, constituting a floor of a quaternary water-bearing layer, in the Vistula valley in the region of KNP occurs from 2 to 54 m a.s.l.The top part of the water-bearing layer has a sandy and sandy-gravel character; the bottom is created by sandy-silt sediments, in places changing to sandy clay and dusty clay [23][24][25].The hydrogeological conditions of KNP have been the subject of several research projects, summarized in several publications [4,5,26,27].Since its resources are significant, this is the main usable aquifer and the major groundwater basin No. 222 [28].The regional groundwater flow occurs to the north, towards the Vistula River.In some regions, there is also a local groundwater circulation, where recharge areas are represented by dunes of higher elevations, whereas swamp areas with their minor watercourses and canals are the drainage base [29].The study area comprises characteristic longitudinal zones of similar geological structure and hydrogeological conditions [4]: the Vistula flood plain terrace, northern and southern marsh zones, northern and southern dune zones, and the accumulative-erosional Warsaw-Błonie terrace (Figure 1).

Analysis of Changes in the Depth to GroundwaterTable
The Vistula River Valley is a typical shallow groundwater area.The groundwater levels are under regular monitoring within a stationary monitoring network.Groundwater level measurements have been conducted since 1999 on a fortnight basis in 56 piezometers in KPN.An analysis of changes of the groundwater table level has provided data on the quantity of groundwater resources.They include the determination of: changes over a long period [4,25], seasonal changes of resources [30], changes in the amount of resources during the vegetation period and the periods and sites of possible moisture deficit in the GDE areas [22], and duration of and depth to groundwater table during hydrogeological drought [17].

Analysis of Changes in the Depth to GroundwaterTable
The Vistula River Valley is a typical shallow groundwater area.The groundwater levels are under regular monitoring within a stationary monitoring network.Groundwater level measurements have been conducted since 1999 on a fortnight basis in 56 piezometers in KPN.An analysis of changes of the groundwater table level has provided data on the quantity of groundwater resources.They include the determination of: changes over a long period [4,25], seasonal changes of resources [30], changes in the amount of resources during the vegetation period and the periods and sites of possible moisture deficit in the GDE areas [22], and duration of and depth to groundwater table during hydrogeological drought [17].
Statistical analysis of groundwater levels in Vistula valley included basic statistics for the 1999-2013 period and particular years of period were analyzed and the scope of groundwater levels was defined.The statistics allow for determining the range of long-term changes in the depth to groundwater table and seasonal changes, as well as defining the highest and lowest levels in a long period.The average depth to groundwater table in the piezometers of the monitoring network during the period of 1999-2013 was 1.87 m b.g.l., and the range of changes varied from 1.11 to 2.46 m.The individual piezometers record typical seasonal level fluctuations.The years 2003 and 2011 stand out for a distinctly different trend of seasonal variations in the groundwater table depth.This phenomenon was observed in all zones of the KNP-marsh, dune, and floodplain.In 2003, the average depth was 2.11 m b.g.l., and the amplitude was 0.87 m.The highest groundwater levels were recorded in 2011, when the average depth to groundwater table was 1.35 m b.g.l., and the amplitude was 0.56 m (Table 1).
In the period of 1999-2013, in marsh zones, the groundwater table occurred at average depths ranging 0-1 m in the central parts (25% of the total area) to 1-3 m in the peripheries (46%) (Figure 2).The zones where the groundwater table was above the ground surface comprised ca.6.5% of the Statistical analysis of groundwater levels in Vistula valley included basic statistics for the 1999-2013 period and particular years of period were analyzed and the scope of groundwater levels was defined.The statistics allow for determining the range of long-term changes in the depth to groundwater table and seasonal changes, as well as defining the highest and lowest levels in a long period.The average depth to groundwater table in the piezometers of the monitoring network during the period of 1999-2013 was 1.87 m b.g.l., and the range of changes varied from 1.11 to 2.46 m.The individual piezometers record typical seasonal level fluctuations.The years 2003 and 2011 stand out for a distinctly different trend of seasonal variations in the groundwater table depth.This phenomenon was observed in all zones of the KNP-marsh, dune, and floodplain.In 2003, the average depth was 2.11 m b.g.l., and the amplitude was 0.87 m.The highest groundwater levels were recorded in 2011, when the average depth to groundwater table was 1.35 m b.g.l., and the amplitude was 0.56 m (Table 1).In the period of 1999-2013, in marsh zones, the groundwater table occurred at average depths ranging 0-1 m in the central parts (25% of the total area) to 1-3 m in the peripheries (46%) (Figure 2).The zones where the groundwater table was above the ground surface comprised ca.6.5% of the study area.The greatest depth to groundwater table (>5 m) was found in dune zones and at the boundary with the upland area, attaining a maximum value of 13.5 m. study area.The greatest depth to groundwater table (>5 m) was found in dune zones and at the boundary with the upland area, attaining a maximum value of 13.5 m.There are distinct differences between the groundwater level for the long period and the lowest (2003) and highest (2011) levels in the period of 1999-2013.In 2003, the groundwater table in the marsh zones was at depths of 0-1 m (12%) or 1-2 m (49% of the total area), which indicates drying of organic soils whose maximum thickness does not exceed 1.5 m [27].We have also observed that the There are distinct differences between the groundwater level for the long period and the lowest (2003) and highest (2011) levels in the period of 1999-2013.In 2003, the groundwater table in the marsh zones was at depths of 0-1 m (12%) or 1-2 m (49% of the total area), which indicates drying of organic soils whose maximum thickness does not exceed 1.5 m [27].We have also observed that the areas where the groundwater table occurs at depths of >5 m doubled, and the maximum groundwater table depth was determined to be 14.5 m.
In 2011, the groundwater table was ca.0.54 m above the ground surface on average, and locally up to 2 m.This phenomenon was observed in almost the entire area of the marsh zones, i.e., ca.29% of the total study area.The depth intervals of 0-1 and 1-3 m account for 29% of the area each.The calculated maximum depth was 11.89 m.

Assessment of Recharge Value
Hydrodynamic modeling studies in the Vistula valley were performed to provide a quantitative description of the components of groundwater balance within the Quaternary aquifer and create spatial distribution of groundwater depth.The modelling was carried out using the Visual ModFLOW 4.2.software [31].The model was built under the steady state, defined as the average annual hydrodynamic state in the period of 1999-2008 and it concerned variants related to renaturalization of wetland areas [32,33].Lateral boundaries of this hydrogeological system were delimited along the main rivers of the region: in the north the Vistula River, in the west the Bzura River, in the south along the edge of the accumulative-erosional level, and in the east along the border of Warsaw.The lower boundary was set as the top part of poorly permeable Neogene clays.
The relationship between the aquifer system and its environment was projected as boundary conditions: RIVER type-for Bzura, Vistula and small canals; General Head type-along south boundary with groundwater flow into Vistula valley from Warsaw-Blonie terrace side.The discretization step was set to 100 × 100 m.The calibration process consisted of applying various values of infiltration recharge and underground evaporation with simultaneous verification of flow rates in rivers and canals [32].Then, the hydrodynamic model was updated by applying three calculation variants of infiltration: modeling for average groundwater levels from the period of 1999-2013, for the lowest level (in 2003), and for the highest level (in 2011) [30].The correlation statistics obtained at 49 monitoring points (located in the model area) and calculated head were as follows for all variants: standard error of the estimate 0.023-0.044m; root mean squared 0.158-0.308m; normalized RMS 0.953-1.849%.
The rate of infiltration recharge was calculated from the hydrodynamic model for the period of 1999-2013 (average value) and for the lowest and highest groundwater levels (Figure 3).The greatest recharge rate values are found predominantly in dune areas.In marsh zones, which are characterized by the highest evaporimetric potential (hydrogenic soils), due to the shallow depth to groundwater, intensification of evapotranspiration-related losses causes a so-called negative recharge in the annual cycle.An analysis of seasonal changes [22] has shown that the effective recharge in the marsh areas usually occurs in springtime due to thawing.This is why pollutants can migrate with the percolating water, although, in the annual cycle, the observed losses are higher than the recharge rate.This phenomenon was the reason for covering the areas of negative annual balance by the vulnerability analysis.The recharge rate varies spatially within the range of −89-65 mm/year, and in the individual model variants, the differences in the average values are 6-12 mm/year.For the average groundwater level, the average recharge rate in this area was approximately 21 mm/year, with standard deviation of 35 mm/year.In 2003, the evaporation rate from the groundwater table in the marsh areas was higher than the recharge rate, and remained at a similar level to the average (−50-0 mm/year).In the dune areas, the recharge rate was clearly lower, by ca.20 mm/year.The average recharge rate for this period was 15 mm/year, with an excursion up to 25 mm/year.In 2011, the recharge rate was 3 mm/year.In the dune areas, high recharge rates dominated (>60 mm/year), whereas in the marsh areas, low values prevailed (<−50 mm/year), which was related to intense evaporation from the groundwater table occurring at the shallowest depths as compared to the other variants.
model variants, the differences in the average values are 6-12 mm/year.For the average groundwater level, the average recharge rate in this area was approximately 21 mm/year, with standard deviation of 35 mm/year.In 2003, the evaporation rate from the groundwater table in the marsh areas was higher than the recharge rate, and remained at a similar level to the average (−50-0 mm/year).In the dune areas, the recharge rate was clearly lower, by ca.20 mm/year.The average recharge rate for this period was 15 mm/year, with an excursion up to 25 mm/year.In 2011, the recharge rate was 3 mm/year.In the dune areas, high recharge rates dominated (>60 mm/year), whereas in the marsh areas, low values prevailed (<−50 mm/year), which was related to intense evaporation from the groundwater table occurring at the shallowest depths as compared to the other variants.

Assessment of Vulnerability
The parametric method DRASTIC (indexation of parameters) was used to assess the natural vulnerability to contamination, dedicated to porous water-bearing media.The most important assumptions of this method include the advective migration of conservative pollutants from the ground surface through the vadose zone.DRASTIC is a method often used to assess groundwater vulnerability to a wide range of potential contaminants.It is a standard tool for assessing groundwater vulnerability to contamination, and is applied in many countries [2,10,11,[34][35][36][37][38][39][40].
Calculations of the DRASTIC index (IPZ) are conducted according to the following formula [1]: where R-parameter rating (Table 2); W-parameter weight (Table 2); D-depth to groundwater; R-net recharge; A-aquifer media; S-soil media; T-topography; I-impact of vadose zone; C-hydraulic conductivity of the aquifer.An advantage of this method is that the used parameters comprehensively characterize the analyzed multi-aquifer system by providing information about its both typical hydrogeological elements (groundwater table depth, recharge rate, aquifer lithology, hydraulic conductivity) and environmental elements that include terrain topography and the type of soil profile and its parent lithology of the vadose zone.
Because the circulation conditions are very well known in the study area, the analysis was based on various, highly detailed data (Table 2).They comprised the analysis of sections of research and hydrogeological boreholes, geophysical surveys, maps, including those compiled within the framework of the Protection plan of Kampinos National Park [41].
To calculate the index of vulnerability (IPZ index), a high-resolution (10 × 10 m) digital terrain model was also applied.A separate issue was the analysis of spatial variability of hydraulic conductivity, whose values were determined based on grain-size analyses, pumping tests, and direct studies using the PARAMEX [26] and BAT tests [29].For the calculations of groundwater vulnerability to contamination, the original weight values were taken ( [1]; Table 2).
The assessment of groundwater vulnerability to contamination was performed for three variants of the depth to groundwater table and infiltration recharge: variant A-average groundwater levels from the period of 1999-2013, treated as long-term average values; variant B-the lowest long-term levels from 2003; and variant C-the highest long-term levels from 2011.The vulnerability calculations were performed using GIS software.Continuous spatial data (i.e., soil types, aquifer media etc.) were discretized using ArcGIS 10.3 and the Raster calculator tool.

Depth to Groundwater
The depth to groundwater table was determined for the period of 1999-2013 and for the lowest (in 2003) and highest (in 2011) recorded groundwater levels.
In 1999-2013, the largest area was that where the groundwater table was present at 1-3 m depth (more than 260 km 2 ) and where it occurred at shallow depths down to 1 m b.g.l.(more than 170 km 2 ).In 2003, the groundwater table was found at greater depths; the area where the groundwater table occurred at depths of 1-5 m was 430 km 2 , and the area where it remained below a depth of 5 m was greater than 75 km 2 .In 2011, the depth to groundwater table occurred most frequently at depths down to 1 m (more than 320 km 2 ).Relevant parameter rating values were attributed to the individual depth intervals (Table 3).

Net Recharge
The division of the R parameter into 6 ranks was made on the basis of the recharge variability in particular periods of analysis (Table 4).In 1999-2013, the largest area was that where the recharge was in the range of 40-60 mm/year (over 263 km 2 ), and mostly represented dune area.Second area with negative recharge from −50 to 0 mm/year (over 139 km 2 ) occurred in the marsh zones.In 2003, the recharge took lower values with 290 km 2 area represented by range of 20-40 mm/year.Respectively lower values of the negative recharge were a result of a deeper groundwater occurring.In 2011, the recharge was again bimodal with a clearer selection of the highest values, i.e., exceeding 60 mm/year (109 km 2 ) and lowest values, i.e., below −50 mm/year (120 km 2 ).

Aquifer Media
The aquifer lithology is generally inhomogeneous.Slight variations are observed as regards thin till interlayers, weathered tills or alluvial muds in sands.The aquifer lithology is represented by three main sediment types (Table 5): fluvial and glaciofluvial sands (37% of the area), gravels and eolian sands (32%), and fine-grained sands with alluvial muds (26%).Their distribution is related to the extents of the individual generations of the Vistula River terraces, as well as to the near-surface deposits at the shallow groundwater table (Figure 4).

Soil Media
The marsh zones are characterized by the greatest heterogeneity with respect to soil types.In this area, there are peaty soils, black earth soils, as well as alluvial soils composed of the finest grain fractions.The dune areas are covered mainly with podzols and gleyed podzols with poorly

Soil Media
The marsh zones are characterized by the greatest heterogeneity with respect to soil types.In this area, there are peaty soils, black earth soils, as well as alluvial soils composed of the finest grain fractions.The dune areas are covered mainly with podzols and gleyed podzols with poorly developed soil horizons, especially the humic layer.The soil ranking scheme (Table 5) was based on defining their protective ability against the migration of pollutants into groundwater, resulting from their field water capacity [44].

Topography
The area is characterized by low relief gradients, ranging predominantly up to 2.5%.Gradients higher than 3% were found on dune slopes; they account for only 2.2% of the area.The most common areas are those favoring infiltration, where the land gradient does not exceed 1% (78% of the area), with the T parameter rank of 10 (Table 5).

Impact of Vadose Zone
The I parameter has similar classes as the aquifer lithology (Table 5).The shallow groundwater table occurs within the same lithology, causing the parameters A and I to be very similar.Their high mutual correlation can be the basis for concluding on the necessity of modifying the DRASTIC assessment due to repeated information in the input layers (Table 2).

Hydraulic Conductivity
The hydraulic conductivity values range from 13 to 28 m/day in most of the area (more than 265 km 2 ), with an average of ca. 25 m/day in 46% of the area.Slightly higher values, 29-40 m/day, are reported for the Vistula floodplain (26%).Zones of increased water permeability occur in the western, central and north-eastern parts of the area, covering ca.70 km 2 (Table 5).Low values, below 4 m/day, are found locally in the southeastern part of the area, occupying 2.96 km 2 .

Parameter Sensitivity Analysis
To determine the effect of weights of the individual parameters applied in the DRASTIC method, a single parameter sensitivity analysis was performed [3,6,45].This analysis determines, by comparing the theoretical and effective weights, the possibility of optimizing the results depended on the input data and the relations between them.The "effective" weight of a parameter is obtained with regard to the other parameters of DRASTIC assessment.To assess the "effective" weight of each parameter-W was calculated using the following formula: where W is the effective parameter weight; the other symbols are as in Equation (1).

Vulnerability Validation
Validation of vulnerability is carried out through correlation with chemical substances, which can be commonly found in groundwater, and could be treated as an indicator of groundwater pollution [11,46,47].In order to perform the validation of variant assessments of intrinsic vulnerability, a correlation analysis was carried out.Correlation analysis consisted in searching for a simple (linear) dependence between datasets and for dependence determined as a result of graphical interpretation of a point cloud.As an indicator of validation, the averaged chloride concentrations from 44 piezometers were used.The chlorides were determined every quarter in the period of 2006-2009.The validation of vulnerability variant of the average hydrodynamic state for theoretical and effective weights (the period 1999-2013) was made.

Results and Discussion
The DRASTIC index was calculated for three variants of parameter D: variant A-for the average depth to groundwater table from 1999 to 2013; variant B-for the lowest groundwater levels in a long period in 2003; and variant C-for the highest groundwater levels in a long period in 2011.
The groundwater vulnerability to contamination, determined for the average level from the period 1999-2013, is characterized by a high proportion of the area with a moderately high vulnerability class (IPZ = 151-175), accounting for 41% of the study area, i.e., 231.87 km 2 .The moderately high vulnerability occurs mainly in the dune zones (Figure 5).Despite the shallow depth to groundwater table in the marsh zones, the vulnerability is lower by 1-2 points, ranging from medium (IPZ = 126-150) in the margins of the zone, to low (IPZ = 101-125) in its central parts.The low vulnerability class is determined for the Vistula floodplain.The total area classified as having low vulnerability covers 138.80 km 2 .The high vulnerability class (IPZ > 176) is identified in the area of 5.68 km 2 in some patches of the southern zone and in the suburbs of Warsaw, in the eastern part of the study area.
The assessment of groundwater vulnerability to contamination for the average depth and recharge rates in the years 1999-2013 was the basis for the analysis of vulnerability changes caused by the fluctuations in the depth to groundwater table and recharge rates.The vulnerability determined for average groundwater levels from a long-term period has been defined as the base state.Vulnerability maps, constructed based on average input data values in various vulnerability assessment methods and procedures, are the basis for making spatial development and groundwater protection decisions.In the case of hydrogeological systems with unconfined groundwater table and significant changes in its dynamics, it is necessary to analyze vulnerability changes for various variants.
The vulnerability determined for the lowest average levels has been defined as the vulnerability during hydrogeological drought.In the conditions of hydrogeological drought, the vulnerability is distinctly lower throughout the area.The medium vulnerability class occupies the greatest part of the study area, i.e., 216.26 km 2 (338%), mainly in the dune zones (Figure 5).The moderately high vulnerability covers an area of 190.75 km 2 .Despite the shallow depth to groundwater table in the marsh zones, the vulnerability score is lower by 1-2, ranging from medium in the margins of the zone to low in its central parts.The low vulnerability class covers the Vistula River floodplain.The total area classified as having low vulnerability is 138.80 km 2 .The high vulnerability class is determined for the total area of 2.12 km 2 in a few patches of the southern dune zone and in the suburbs of Warsaw, in the eastern part of the study area.
The assessment results were the basis for analyzing the differences in the areas occupied by the individual vulnerability classes determined for the base state and for the conditions of hydrogeological drought.In general, the class distribution pattern has not changed significantly.The most distinct changes refer to the size of highly vulnerable areas (a decrease from 0.99% to 0.37%) and to the increase in the very low class area (from 0.21% to 0.32%).The analysis of spatial variations in the IPZ index values shows that its greatest decrease in relation to the average value was 15, and the average decrease was 3.19.Considering this change with regard to the size of the classes (IPZ = 25), the changes can be deemed insignificant.
In variant C, for the highest groundwater levels in 2011, the vulnerability is close to the base level.The greatest area is covered by the moderately high vulnerability class, which is 221.75 km 2 of the study area, predominantly in the dune zones (Figure 5).Highly vulnerable areas occupy 197.99 km 2 .The low vulnerability class occurs chiefly in the Vistula floodplain, and the total area classified as having low vulnerability is 141.11 km 2 .The high vulnerability class occupies the total area of 9.29 km 2 , in a few patches of the southern dune zone and in the suburbs of Warsaw-in the eastern part of the study area.When analyzing the differences between the base vulnerability and the flood risk vulnerability (Table 6), a slight increase in the area of low vulnerability class and nearly a twice increase in the area of high vulnerability class (from 0.99% to 1.62%; Figure 5) can be noticed.The analysis of spatial differences in the IPZ values has shown an average increase by merely 0.3, with a maximum increase by 14.5, although there are areas where the IPZ decreased (by up to 13).In variant C, for the highest groundwater levels in 2011, the vulnerability is close to the base level.The greatest area is covered by the moderately high vulnerability class, which is 221.75 km 2 of the study area, predominantly in the dune zones (Figure 5).Highly vulnerable areas occupy 197.99 km 2 .The low vulnerability class occurs chiefly in the Vistula floodplain, and the total area classified as having low vulnerability is 141.11 km 2 .The high vulnerability class occupies the total area of 9.29 km 2 , in a few patches of the southern dune zone and in the suburbs of Warsaw-in the eastern part of the study area.When analyzing the differences between the base vulnerability and the flood risk vulnerability (Table 6), a slight increase in the area of low vulnerability class and nearly a twice increase in the area of high vulnerability class (from 0.99% to 1.62%; Figure 5) can be noticed.The A single parameter sensitivity analysis was carried out and, after determination of new weights, the IPZ index was recalculated for the above-mentioned three variants dependent on the groundwater levels and recharge rates (Figures 2 and 3).The greatest increase in the effective weight in comparison with the theoretical weight refers to parameter D; the effective value exceeds the theoretical one by more than 50% (Table 7; Figure 6).Similar values of effective weights were also obtained in the sensitivity analysis concerning Vistula flood plain itself [48], which means that this dependence is characteristic for the whole area covered by the analysis.For parameter I, which has the theoretical weight of 5, like parameter D, the calculated effective weight is lower and determined at 4.4.For parameter R, which is characterized by the theoretical weight of 4, the effective weight is much lower and it was calculated on 2.4.Lower effective weights have also been obtained for parameter A, from 3 to 2.7, and for parameter C, from 3 to 2.3 (Table 7).The calculations made by the single parameter sensitivity analysis confirm the role of the depth to groundwater table in the assessment of groundwater vulnerability to contamination using the DRASTIC method.In relation to the vulnerability classes determined earlier, there was an increase in the IPZ index in all cases (Figure 7; Table 8).The percentage of the area of the lowest vulnerability classes, in particular the low class, decreased by ca.23-24%.The area occupied by the moderately high and very high classes increased by several percent.An increase in the very high class area is lower (by 6.1%) only for the low groundwater level variant (in 2003).Moderately high and very high classes dominate in the dune areas, whereas the axial parts of the marsh areas are covered by the medium class (locally passing into the low class), with the moderately high class at the margins.The largest area of the low and very low classes occurs in the south, along the margin of the Błonie Level.The change in the DRASTIC index value for the average level during the period was from −1 to 26, at the average value indicating an increase, after optimization, by 14.4.In 2003, the index change varied from 2 to 27, with an average of 13.9, while in 2011, the range of changes was 0-27, with an average of 16.The greatest increase in the index value was found for the marsh areas, which means that the contribution of parameter D is underestimated when the theoretical weights remain at the same level, and the obtained vulnerability classes are lower than those inferred from the actual effect of the depth on the ultimate vulnerability assessment.No changes are observed in the dune areas, where the vulnerability change is low (Figure 7).The study area is characterized by low concentrations of chlorides, which reach values from 4.2 to 111 mg/L, with an average of 28.2 mg/L (median: 17.6 mg/L).The lowest concentrations occur on dune areas, while the highest are in drainage zones near the canals.Validation of vulnerability as an expected high dependence of chloride concentration and IPZ index was found as the reverse.For low values of Cl − ions, the highest IPZ indexes were calculated (Figure 8).This results from the lack of source of groundwater contamination and intensity of infiltration.Intense recharge in the area of   The study area is characterized by low concentrations of chlorides, which reach values from 4.2 to 111 mg/L, with an average of 28.2 mg/L (median: 17.6 mg/L).The lowest concentrations occur on dune areas, while the highest are in drainage zones near the canals.Validation of vulnerability as an expected high dependence of chloride concentration and IPZ index was found as the reverse.For low values of Cl − ions, the highest IPZ indexes were calculated (Figure 8).This results from the lack of source of groundwater contamination and intensity of infiltration.Intense recharge in the area of dune belts causes a slight transformation of precipitation water.In drainage zones, due to the long residence time of water in the system, processes such as evaporation result in much higher chloride concentrations.dune belts causes a slight transformation of precipitation water.In drainage zones, due to the long residence time of water in the system, processes such as evaporation result in much higher chloride concentrations.
The results of the correlation analysis indicate a better adjustment of the trend line to the results of the DRASTIC method with the original weights, where lesser importance is attributed to depth, and higher to the recharge and lithology of the aquifer and the vadose zone.The results of the correlation analysis indicate a better adjustment of the trend line to the results of the DRASTIC method with the original weights, where lesser importance is attributed to depth, and higher to the recharge and lithology of the aquifer and the vadose zone.

1.
The groundwater vulnerability assessment in the Vistula Valley was performed for three variants: base variant A-average depth to groundwater from the period of 1999-2013, hydrogeological drought variant B (lowest long-term depth-2003), and flood risk variant C (greatest long-term depth-2011).In variant A, the moderately high and medium vulnerability classes occupy the greatest area.During hydrogeological drought, the largest area is occupied by the medium vulnerability class.During the potential flood risk in the river valleys, the greatest area is occupied by the moderately high vulnerability class.All these variants reveal a spatial variability in the distribution of the individual classes.

2.
Based on sensitivity analysis, depth to water table was the most effective parameter responsible for highest variation in the vulnerability index.The effective weight of parameter D is much higher than the theoretical value.In all variants, the calculated IPZ index indicates (after optimization of the parameter values) the greatest increase in the area of the moderately high and high classes, while the area of the low class decreased.The greatest changes in relations to the theoretical values occurred in variant B for the lowest groundwater level variant.

3.
The studies confirm that the assessment of groundwater vulnerability to contamination in areas of shallow depth to groundwater table, which are GDE, should be made for a number of variants, with a particular focus on the highest and lowest observed groundwater levels.Groundwater vulnerability changes are also associated with sensitivity analysis and calculated effective weights.

4.
This study suggests that DRASTIC is an effective tool for the groundwater vulnerability assessment studies and can be used for the prioritization of the susceptible areas to avoid contamination as well as a frequent and detailed monitoring of pollution potential in the high and medium vulnerability zones.There is also a need for proper planning in areas of shallow groundwater and considerable changes.Vulnerability maps for various variants of input data can be the basis for plans of land management and flood protection, as well as for prediction of climate change effects.

Figure 1 .
Figure 1.3D view of study area location (vertical exaggeration 5.0; view from west; gridlines every 2 km).

Figure 1 .
Figure 1.3D view of study area location (vertical exaggeration 5.0; view from west; gridlines every 2 km).

Figure 8 .
Figure 8. Correlation of chloride concentration and IPZ index in selected piezometers.

Table 1 .
Main statistics of groundwater depth based on monitoring observations in all the piezometers and in selected piezometers of dune (P11), marsh (P12), and Vistula flood plain (P14).

Table 2 .
DRASTIC parameters and sources of data used.

Table 3 .
Groundwater depth and its range in selected years.

Table 4 .
Net recharge and its range in selected years.

Table 5 .
Characteristics and range of A, S, T, I, C parameters.

Table 6 .
Statistics of IPZ index and occurrence of vulnerability classes on research area.

Table 7 .
Original and effective weight of parameters.
Comparison of original and effective weights of parameters.

Table 8 .
IPZ changes after optimization by single parameter sensitivity analysis.

Table 8 .
IPZ changes after optimization by single parameter sensitivity analysis.