Evaluation of the Inﬂuence of Farming Practices and Land Use on Groundwater Resources in a Coastal Multi-Aquifer System in Puck Region (Northern Poland)

: This study focuses on the modeling of groundwater ﬂow and nitrate transport in a multi- aquifer hydrosystem in northern Poland, adjacent to Puck Bay (Baltic sea). The main goal was to investigate how changes in land use and farming practices may a ﬀ ect groundwater recharge and submarine groundwater discharge (SGD) to the sea and the associated N-NO 3 ﬂuxes. An integrated modelling approach has been developed, which couples the SWAT hydrologic model, MODFLOW-NWT groundwater ﬂow model, and MT3DMS transport model. Transient simulations were performed for a 10 y period, assuming 10 di ﬀ erent scenarios of land use (farming, grassland, forest) and crop types. Both recharge and SGD showed a distinct pattern of seasonal time variability. In terms of the average ﬂow rates, the e ﬀ ect of varying crop type was more signiﬁcant than that of land use change, with the minimum recharge and SGD corresponding to winter wheat and the maximum for peas and potatoes. Nitrate loads were strongly a ﬀ ected by both land use and crop type, with minimum values obtained for grassland and maximum values for canola.


Introduction
The development of integrated surface-subsurface flow and transport models is often based on the coupling of existing computer codes, capable of simulating processes in various hydrological compartments. The combination of SWAT hydrologic model [1] with MODFLOW/MT3D groundwater flow and transport models [2,3] has been used by a number of authors over the last 20 years . Such an approach enhances the capacities of both models because the groundwater component of SWAT is simplified and does not offer the accuracy required in many applications, especially involving complex hydrological conditions. On the other hand, SWAT is well-suited for computations of recharge rates and the associated nitrate loads varying in space and time over the watershed. Coupling between computer codes can be implemented as one-way sequential coupling, where spatially distributed recharge rates and nitrate loads from SWAT (and possibly also river stages) are passed to MODFLOW and MT3DMS, which calculate groundwater flow and transport independently from the SWAT groundwater subroutine. Alternatively, the groundwater procedure in SWAT can be replaced by the MODFLOW code, in principle allowing for a more consistent representation of surface-subsurface water interaction, at the cost of increased complexity of the computer code. Recently, a fully coupled SWAT-MODFLOW-RT3D software was developed in the form of a single executable program [15].
In previous works SWAT was coupled with MODFLOW for watersheds varying in size from a few km 2 [24] to several thousand km 2 [19]. However, there are only a few studies that used the combination of SWAT, MODFLOW, and MT3D to evaluate the quantity and quality of the submarine groundwater discharge (SGD) from land aquifers to coastal marine waters. SGD has received increasing attention of scientists due to its role played in water and nutrient cycles [25][26][27][28][29]. Galbiati et al. [7] developed an integrated model consisting of the SWAT, MODFLOW, MT3D, and stream water quality model QUAL2E in order to estimate nutrient load discharged to an Adriatic coastal lagoon. In their case, however, the nutrient influx was caused mainly by pumping from channels, not SGD, because the land surface was mostly below the sea level. More recently, [24] applied recharge data from SWAT as an input to the MODFLOW/MT3D model describing local groundwater flow in a small coastal watershed on Samoa. On the other hand, there is only limited experience in applying SWAT to complex multi-aquifer groundwater models [14,19] and, to the best of our knowledge, there were, until now, no attempts to use SWAT coupled with a multi-aquifer groundwater model to estimate SGD. The research reported here has been motivated by the need to predict nutrient loads in submarine groundwater discharge to the Puck Bay (northern Poland), which is a shallow part of the Baltic Sea, very vulnerable to eutrophication ( Figure 1) [28,[30][31][32]. Our study is a part of a larger project, which combines SWAT, MODFLOW, and MT3D with a hydrodynamic and ecological model of the Puck Bay (EcoPuckBay [33,34]). The integrated set of models can be used as a decision supporting tool and constitutes a part of the web-based service WaterPUCK (waterpuck.pl), which gives local stakeholders (administration and farmers) tools for predicting changes in the quality and quantity of local water resources. Below, we list the objectives of our work and show how they address the existing gaps of knowledge: 1.
Development of a numerical model describing transient groundwater flow and nitrate transport in the multi-aquifer system adjacent to the Puck Bay. Earlier groundwater modeling studies in this area focused only on the rate of SGD, without contaminant transport [35][36][37].

2.
Coupling of groundwater model with SWAT hydrological model to obtain detailed information on groundwater recharge rate and nitrate loads leached from soil. SWAT has already been used in a number of studies on the hydrological balance and nutrient transport in the Baltic region [38][39][40][41][42][43], including the Puck Bay region [44,45]. However, there has been no attempt to apply SWAT-MODFLOW coupling on the Polish Baltic coast.

3.
Evaluation of the time variability of the groundwater recharge rate and SGD rate and the associated nitrate loads under different land use scenarios. While the impact of crop type and land use on water and nutrient fluxes has been extensively studied for different watersheds using SWAT, there have been no such investigations in the context of submarine groundwater discharge. In this work, we evaluate 10 scenarios, corresponding to six types of crops and four types of land use.

Study Area
The area under investigation is located in northern Poland ( Figure 1). The land represents a typical young glacial landscape with relatively high relief. It consists of isolated fragments of a moraine plateau separated from each other with deeply cut ice marginal valleys. The entire area is covered with Quaternary deposits. The thickness of the sediments ranges from 40 m in the east to 95-100 m in the west. Such a significant difference is probably a result of the varied land surface or large range in altitude of the top of older layers [46]. The Quaternary deposits consist of moraine glacial till with sand and gravel layers, and glacifluvial or river sand and silty sand in the valleys. The dominant soil types are sandy loam (glacial till), sandy loam covered by loamy sand (weathered glacial till), sand (of glacifluvial origin), and peat (in larger river valleys).

Study Area
The area under investigation is located in northern Poland ( Figure 1). The land represents a typical young glacial landscape with relatively high relief. It consists of isolated fragments of a moraine plateau separated from each other with deeply cut ice marginal valleys. The entire area is covered with Quaternary deposits. The thickness of the sediments ranges from 40 m in the east to 95-100 m in the west. Such a significant difference is probably a result of the varied land surface or large range in altitude of the top of older layers [46]. The Quaternary deposits consist of moraine glacial till with sand and gravel layers, and glacifluvial or river sand and silty sand in the valleys. The dominant soil types are sandy loam (glacial till), sandy loam covered by loamy sand (weathered glacial till), sand (of glacifluvial origin), and peat (in larger river valleys).
Groundwater forms a complex system with several aquifer units ( Figure 2). Shallow groundwater occurs in form of small perched aquifers and in sand lenses enclosed in moraine deposits (Qp in Figure 1). These aquifers are situated at depths up to 5-10 m below the ground level. They do not play an important role in water supply, but they are exploited locally by farms and summer houses. Two deeper Quaternary aquifers span most of the area: The upper (sub-moraine) aquifer (Q1) and the lower inter-moraine aquifer (Q2). They were formed in glacifluvial deposits (sand and gravel) separated with a layer of moraine till. These two large aquifers are hydraulically connected. On the major part of the area, they are confined (piezometric head reaches 45 m above mean sea level), except areas in the valleys where the till cover was eroded and the upper aquifer is exposed. In those places, the upper aquifer is unconfined and indirectly connected with surface waters. Shallow groundwater is recharged by the infiltration of rainwater. Deeper aquifers are influenced by inflow from the west and take the seepage water from the layers above. The entire hydrosystem is drained mainly by the Baltic Sea (Puck Bay), either directly via SGD or indirectly via streams and rivers. Aquifers Q1 and Q2 are the main source of water supply in the region. Groundwater forms a complex system with several aquifer units ( Figure 2). Shallow groundwater occurs in form of small perched aquifers and in sand lenses enclosed in moraine deposits (Qp in Figure 1). These aquifers are situated at depths up to 5-10 m below the ground level. They do not play an important role in water supply, but they are exploited locally by farms and summer houses. Two deeper Quaternary aquifers span most of the area: The upper (sub-moraine) aquifer (Q1) and the lower inter-moraine aquifer (Q2). They were formed in glacifluvial deposits (sand and gravel) separated with a layer of moraine till. These two large aquifers are hydraulically connected. On the major part of the area, they are confined (piezometric head reaches 45 m above mean sea level), except areas in the valleys where the till cover was eroded and the upper aquifer is exposed. In those places, the upper aquifer is unconfined and indirectly connected with surface waters. Shallow groundwater is recharged by the infiltration of rainwater. Deeper aquifers are influenced by inflow from the west and take the seepage water from the layers above. The entire hydrosystem is drained mainly by the Baltic Sea (Puck Bay), either directly via SGD or indirectly via streams and rivers. Aquifers Q1 and Q2 are the main source of water supply in the region.
Proximity to the Baltic Sea results in a specific marine climate of this region typically with moderate winters and mild summers. The annual temperature is about 7.4 °C and the average precipitation in the period of 2000-2009 was 620 mm/y.
The population of the area is about 25,000 people, of which 11,000 live in the largest town, Puck. Land use is divided between agricultural (51%), forests (28%), permanent meadows and grasslands (12%, mostly on peat soils in river valleys), and others, mostly low-density urban (9%). The main crops are (in order of decreasing importance): Winter wheat, winter triticale (cross-over between Proximity to the Baltic Sea results in a specific marine climate of this region typically with moderate winters and mild summers. The annual temperature is about 7.4 • C and the average precipitation in the period of 2000-2009 was 620 mm/y. The population of the area is about 25,000 people, of which 11,000 live in the largest town, Puck. Land use is divided between agricultural (51%), forests (28%), permanent meadows and grasslands (12%, mostly on peat soils in river valleys), and others, mostly low-density urban (9%). The main crops are (in order of decreasing importance): Winter wheat, winter triticale (cross-over between wheat and rye, typical for Poland), silage corn, canola, mixed spring cereals, potatoes, and peas (Pisum). Most of the farmers use crop rotation [47].

SWAT Model
A single SWAT model has been built for the investigated area, comprising the watersheds of 3 streams flowing to Puck Bay (Gizdepka, Płutnica, Błądzikowski Potok), as well as a small part of the watershed of Reda river, adjacent to its outflow to the bay ( Figure 1). The model has been set-up using the QGIS SWAT interface. The area was divided into 17 subbasins and 353 hydrological response units (HRU-basic units for SWAT calculations). HRUs were defined based on the land use map and soil Water 2020, 12, 1042 5 of 17 map. Screening tests showed that including ground slope as a criterion significantly increases the number of HRUs, without affecting the results of simulations. Weather data for the period 2000-2009 were obtained from the station inŻelistrzewo, in the SE part of the area (daily precipitation) and from Climate Forecast System Reanalysis [globalweather.tamu.edu] (daily values of humidity, wind, temperature, and radiation). Potential evapotranspiration was estimated with the Penman-Monteith model. Soil parameters (available water content, vertical hydraulic conductivity, etc.) were taken as average values in each textural class, as reported in the Polish literature [48]. The agricultural management practices for crops and grassland were based on the results of a poll carried out among local farmers [47]. In all agricultural areas, the same pattern of crop rotation was specified, including silage corn, canola, and winter wheat. Similarly, all grassland was assigned the same set of practices (permanent plant cover, Fescue species, hay cut twice a year). Forests were modeled as permanent plant cover. Pine was selected as the representative species.
Due to the lack of sufficient data on stream discharges, the calibration of SWAT has been performed by a manual trial-and-error procedure, using the available data on the components of the hydrological budget and crop yield (soft calibration). The SWAT model was calibrated in parallel with the steady-state groundwater flow model, which provided an additional constraint in the calibration process. As the first objective for SWAT calibration, we set the limits of the calculated groundwater recharge in any HRU from 3% to 30% of precipitation, which is commonly accepted for the climate and soil conditions in Poland [49,50]. The lower limit corresponded to HRU's on peat and heavy clay soils, while the upper limit -to clean glaci-fluvial sands. As the second objective, we set the value of the average yearly evapotranspiration, which was reported to be in the range from 450 to 495 mm for a location in Poland with very similar climate and land use conditions (629 mm precipitation, canola with different fertilization schemes) [51]. The third objective was the annual biomass production in the forests, which should be close to 7 metric tons of dry mass per hectare per year [54,55].
In the calibration process, the following parameters were adjusted: SCN curve numbers (increased), maximum canopy storage (CANMX, increased in forests), and REVAP (increased in forests), as well as the parameters describing pine growth in the "plant.dat" file. A satisfactory agreement was obtained between the reference values (calibration goals) and the results obtained from the calibrated model (Table 1). In the same table, we report additional validation checks, performed on model outputs that were not used previously in the calibration procedure. The yield of basic crop types (winter wheat, canola, and silage corn) and the biomass harvested as hay on grasslands were compared to the data obtained from local farmers, from the statistics for the Pomeranian region in Poland and reported in the literature (all values given in Table 1 represent dry mass). The model predictions are quite accurate for canola, while there is some underestimation of yield for silage corn and some overestimation for winter wheat. Moreover, both the total runoff and the surface runoff/total runoff ratio are in agreement with previous studies on the hydrology of the region. Overall, the performance of the model versus the available data can be considered satisfactory. In groundwater simulations described below, we used, as input data, the values of groundwater recharge and N-NO 3 load reaching the groundwater table, obtained from SWAT simulations for the period 2000-2009, with a 2 y warm-up period (1998)(1999). Note that SWAT uses a simplified procedure to describe the movement of water and nitrate through the vadose zone between the bottom of the soil profile and the groundwater table [1]. The resulting estimations of the vadose zone time lag may be inaccurate, but at the current stage of model development, this simplification was accepted.

Groundwater Flow Model
Groundwater flow simulations were carried out with MODFLOW-NWT numerical code [57]. The model was set-up with GMS 10.4 [58] and ModelMuse [59] graphical interfaces, while the FloPy [60] library was used to prepare data for transient simulations.
In line with the finite difference discretization of MODFLOW, we used a rectangular grid with 6 layers, 410 rows, and 296 columns. The total number of cells was 728,160, of which 340,496 were active. The spacing in horizontal directions was uniform and equal to 50 m in order to capture the topographic features and minimize discretization errors. The 6 layers were defined as follows: (I) Sandy, perched aquifers, occurring locally (Qp), (II) low-permeability till deposits, (III) sub-moraine aquifer in glacifluvial sand (Q1), (IV) low permeable till layer, (V) inter-moraine aquifer in glacifluvial sand (Q2), (VI) low-permeability layer representing Quarternary glacial till and Neogene clay and silt. The range of hydraulic conductivity values for each layer is presented in Table 2.
In all steady and transient simulations, we used the same set of boundary conditions and source terms, assumed to be constant in time, except for the recharge rate. The sea boundary was represented by a constant head boundary in the Q1 aquifer (direct hydraulic contact with water in Puck Bay) and by a general head boundary (GHB, third type) in the Q2 aquifer, which represented the outflow through a distant outcrop at the sea bottom (the discharge zones are situated in the bottom of the Puck Bay, at a distance of 2-3 km from the coast). The GHB boundary was also used to simulate lateral flow across the land part of the boundary in Q1, Q2, and those perched aquifers (Qp), which were intersected by the model boundary. Wells were assumed to be constantly in operation and the pumping rates were identified based on the available reports from groundwater administration. Rivers were represented by the third-type boundary conditions (RIVER package). The river stages were assumed to be constant in time and represented average values from measurements and earlier reports. The results obtained from SWAT were processed by scripts written in the Python language using the FloPy library, in order to set the recharge values in the MODFLOW model. In steady-state simulations, we used the average values in each HRU from the period of 2000-2009. For transient simulations, monthly averages were calculated for each HRU from SWAT daily results. We followed the approach described by [15] to transfer HRU-based SWAT results to the grid-based MODFLOW input.
Calibration of the groundwater flow model was based on the steady-state solution, representing average conditions in the period of 2000-2009. We used groundwater head measurements from 95 points in the area. A good agreement was obtained (mean error 0.29 m, mean absolute error 1.10 m, and standard error 1.79 m). The comparison of calculated and measured groundwater head values for all three aquifers (Qp, Q1, and Q2) is shown in Figure 3. In order to validate the model, a number of simulations based on historical records of pumping tests were performed. They confirmed that the model was calibrated adequately. Calibration of the groundwater flow model was based on the steady-state solution, representing average conditions in the period of 2000-2009. We used groundwater head measurements from 95 points in the area. A good agreement was obtained (mean error 0.29 m, mean absolute error 1.10 m, and standard error 1.79 m). The comparison of calculated and measured groundwater head values for all three aquifers (Qp, Q1, and Q2) is shown in Figure 3. In order to validate the model, a number of simulations based on historical records of pumping tests were performed. They confirmed that the model was calibrated adequately.

Nitrate Transport Model
Based on the groundwater flow model developed in MODFLOW, a simulation of N-NO3 transport was carried out. MT3DMS numerical code was applied to solve the advection-reactiondispersion equation. Transient calculations were performed with the third-order total-variational diminishing (TVD) numerical method for the advective term, while for steady-state solution, the standard finite difference discretization was applied.
Porosity and longitudinal dispersivity were assigned depending on the soil material. The values were selected based on the literature. Longitudinal dispersivity was assumed to be 1.5 m in permeable layers [61] and 0.75 m in the case of low-permeability layers [62]. Values of porosity were assigned equal to 0.2 in sand and gravel and 0.1 in glacial till. Denitrification of nitrates was described with the first-order kinetic reaction [63][64][65][66]. We used a constant reaction rate of biodegradation equal to 1x10 -5 1/h, which was established by model calibration.
In the model, we considered only diffuse sources of nitrates. Concentrations of N-NO3 in water reaching the groundwater table were calculated from the results of SWAT simulation for the years of 2000-2009 for each HRU. As in the case of groundwater recharge, the transfer of HRU-based SWAT data to grid-based MT3DMS data was carried out with the aid of Python scripts, making use of the FloPy library. In the areas of the groundwater model not covered by SWAT, we used a concentration value representing an average taken over the whole SWAT model area in the considered time period.

Nitrate Transport Model
Based on the groundwater flow model developed in MODFLOW, a simulation of N-NO 3 transport was carried out. MT3DMS numerical code was applied to solve the advection-reaction-dispersion equation. Transient calculations were performed with the third-order total-variational diminishing (TVD) numerical method for the advective term, while for steady-state solution, the standard finite difference discretization was applied.
Porosity and longitudinal dispersivity were assigned depending on the soil material. The values were selected based on the literature. Longitudinal dispersivity was assumed to be 1.5 m in permeable layers [61] and 0.75 m in the case of low-permeability layers [62]. Values of porosity were assigned equal to 0.2 in sand and gravel and 0.1 in glacial till. Denitrification of nitrates was described with the first-order kinetic reaction [63][64][65][66]. We used a constant reaction rate of biodegradation equal to 1 × 10 −5 1/h, which was established by model calibration.
In the model, we considered only diffuse sources of nitrates. Concentrations of N-NO 3 in water reaching the groundwater table were calculated from the results of SWAT simulation for the years of 2000-2009 for each HRU. As in the case of groundwater recharge, the transfer of HRU-based SWAT data to grid-based MT3DMS data was carried out with the aid of Python scripts, making use of the FloPy library. In the areas of the groundwater model not covered by SWAT, we used a concentration value representing an average taken over the whole SWAT model area in the considered time period.
The model of nitrate transport was calibrated manually by adjusting the reaction rate in order to obtain a satisfactory agreement between the calculated N-NO 3 concentrations and the results of chemical analysis of the samples taken from the aquifers. Calibration was based on a steady-state flow and transport simulation, using average values of recharge and nitrate loads from the period 2000-2009 and the current land use and agricultural practices. The range of N-NO 3 concentrations in layer Q1 obtained from the numerical simulation was 3.8 × 10 −6 -62.88 mg/dm 3 , while the range of measured concentrations was <0.02-20.78 mg/dm 3 (0.02 is the detection limit). The values obtained from numerical modeling for layer Q2 were from 0.0012 to 1 mg/dm 3 , while the measured values ranged between <0.02 and 0.23 mg/dm 3 .

Scenarios for Transient Flow Simulations
In order to investigate the effects of different land use on the time variability of groundwater recharge and groundwater discharge to the Puck Bay, and the associated N-NO 3 loads, we defined 10 scenarios (see Table 3). Scenario S1 (baseline) corresponds to the current land use and agricultural practices. In agricultural areas, crop rotation is used with alternating growing of winter wheat, silage corn, and canola, which are the main crops in the area.
In scenarios S2-S7, we assumed that the land use structure corresponds to the current state, but there is only one type of crop grown on agricultural lands, without any rotation. The crop types are as follows: S2: Winter wheat (also representing winter triticale), S3: Silage corn, S4: Canola, S5: Mixture of spring cereals (represented by barley), S6: Potatoes, S7: Peas (Pisum).
Scenarios S8-S10 represent major shifts in land use structure: S8: Whole area covered by agriculture, S9: Whole area covered by meadows, S10: Whole area covered by forest (in each case, we kept the current location of urban areas). While no such large changes are predicted in the near future, these scenarios can be considered as limit cases, defining the maximum range of expected variations in groundwater fluxes and nitrate loads.

Steady-State Flow Simulation
The components of the steady-state groundwater budget are shown in Figure 2. According to the model results, the main components of the budget are recharge from precipitation, external lateral inflow to Q1 and Q2, exchange with streams, and submarine groundwater discharge. The considered groundwater system is typical for young glacial areas, with flow taking part in a multiple stacked aquifer, partly connected with each other. A typical feature of such systems is significant vertical seepage through glacial till layers, separating the aquifers. Even though glacial till is often considered a weakly permeable deposit, it plays an important role in the circulation of groundwater on glacial areas. These results are consistent with previous research, e.g., [50].
Submarine groundwater discharge constitutes a major component in groundwater budget. The amount of water discharged from the lower aquifer Q2 is approximately twice as large as in the case of the upper aquifer Q1. This can be explained by significant lateral inflow received by Q2 aquifer from the western boundary. In our model, the total groundwater discharge to the sea from Q1 and Q2 is 1286 m 3 /h, which corresponds to 92 m 3 /h per kilometer of the coastline. The flow rate varies from 58 m 3 /h/km in the upland areas to 152 m 3 /h/km in the river valleys. These values are in a good agreement with the results of previous studies related to groundwater outflow to the Puck Bay [35][36][37]67]. Ref. [35] reported the SGD rate in the same area ranging from 21 to 165 m 3 /h/km, with an average of 98 m 3 /h/km, while [37] estimated the SGD rate in the Reda river valley as 165 m 3 /h/km.

Transient Simulations of Flow and Nitrogen Transport
A summary of the main results obtained from transient flow simulations is presented in Table 4. We report the values of groundwater discharge and the corresponding load of N-NO 3 leached from soil to groundwater, as well as the rate of groundwater discharge to Puck Bay and the associated loads of N-NO 3 , separately from Q1 and Q2 aquifers. All values represent averages for the whole simulation period (120 months).
Comparison of the average values of recharge and SGD for different scenarios indicates that the effect of changing the crop type within the current agricultural area (S2-S7) can be more pronounced than large-scale changes in land use (S8-S10). SWAT simulations predicted the minimum recharge for winter wheat (S2: 62 mm/y) and the maximum for peas (S7: 92 mm/y), which corresponds to a relative difference of about 50%. In contrast, the difference between land use types is smaller, with the maximum value of 85 mm/y corresponding to grassland (S9) and the minimum value of 71 mm/y to the forest (S10). The average SGD flow rates for various scenarios show a similar variability pattern as the values of recharge. The relative difference between maximum value (S7) and minimum value (S2) is almost 14% in the Q1 aquifer and 11% in the Q2 aquifer.
The yearly nitrate loads obtained from SWAT for various crops (Table 3) are more varied than recharge values, ranging from 7.4 kg/ha/y for winter wheat to 30.5 kg/ha/y for canola. While in the case of winter wheat, the minimum load corresponds to minimum recharge, the opposite is not true for canola. Therefore, there is no simple relationship between the amount of recharge and nitrate load to groundwater, because the latter one strongly depends on crop type and agricultural practices. Overall, the amount of N leached to groundwater according to SWAT simulations is similar to the values reported in the literature for Poland [68] (5-15 kgN/ha/y for grasslands, 5-10 for forests, 20-45 for crops).
Considering nitrate loads entering Puck Bay from SGD, one can see that the relative magnitude of the load for various scenarios corresponds to the relative magnitude of nitrate leaching from soil. Variation in the N-NO 3 load carried by SGD between scenarios is much larger that the variation in the SGD flow rate. The major part of the load (about 90%) is discharged from the upper aquifer Q1, despite the fact that its flow rate is about two times smaller than in Q2. This can be easily explained, as Q1 is closer to the ground surface and directly receives NO 3 leached from soil. The results of transient simulations show significant time variation of the components of the groundwater budget and the associated nitrate loads. In Figure 4, Figure 5, and Figure 6, we show average monthly water and nitrogen fluxes in recharge and SGD, for scenarios S1 (baseline), S2 (winter wheat: Low recharge and nitrate load), and S3 (silage corn: High recharge and nitrate load). In most cases, a consistent pattern of seasonal changes can be distinguished, with maximum values in late winter/early spring and minimum values in early autumn. Such a variability is typical for Polish climatic conditions-the maximum corresponds to the period of thawing snow and the minimum to the late vegetation season. An exception can be observed in the very dry years of 2005 and 2006, for which there is no early spring maximum in 2006 in S1 and S2, due to higher water demands of plants in these scenarios.
Seasonal changes in SGD (Figures 4-6) correspond approximately to seasonal changes in groundwater recharge; however, the relative differences in monthly values are significantly smaller.
This can be seen particularly in the discharge from Q2, which is less affected by variations in recharge and is driven to a larger extent by lateral inflow from the western boundary. The relative fluctuations of SGD from Q1 reach 120% of the lowest value, while in Q2, they do not exceed 31% of the lowest value.
Concentrations of N-NO 3 in SGD flux obtained from numerical simulations were compared to the results of chemical analysis of groundwater samples (Table 4). In both cases, the values are similar considering that the SGD concentrations from numerical modeling represent the entire area of investigations, and field research was conducted at points distributed across the site. Seasonal changes in SGD (Figures 4-6) correspond approximately to seasonal changes in groundwater recharge; however, the relative differences in monthly values are significantly smaller. This can be seen particularly in the discharge from Q2, which is less affected by variations in recharge and is driven to a larger extent by lateral inflow from the western boundary. The relative fluctuations of SGD from Q1 reach 120% of the lowest value, while in Q2, they do not exceed 31% of the lowest value.
Concentrations of N-NO3 in SGD flux obtained from numerical simulations were compared to the results of chemical analysis of groundwater samples (Table 4). In both cases, the values are similar considering that the SGD concentrations from numerical modeling represent the entire area of investigations, and field research was conducted at points distributed across the site.    As shown in scenarios S8-S10, nitrate loads leaching from soil and transported by SGD are much more sensitive to land use changes than recharge and SGD flow rates. Our simulations show that converting forests and meadows to agriculture (S8) results in nitrate loads increased by about 75% with respect to the baseline scenario. This is consistent with the increase in agricultural area from 51% to 91% of the watershed. The large nitrate load is due to growing corn and canola, as part of crop rotation, assumed in the scenario.
Conversely, if crops and forests are replaced by grassland (S9), the nitrate loads are reduced by about six times compared to the baseline scenario. Grasslands are characterized by a high uptake of nitrogen, promoted by cutting hay twice a year, and the fertilized doses are also smaller.
The introduction of forests to the whole area (S10) also leads to the reduction of nitrate loads, albeit much less significant than in scenario S9. This result is caused by a smaller uptake of nitrogen in the grown-up forest as compared to growing grass. Even though there is no fertilization in the forest, a significant amount of nitrogen is released from peat soils. In other scenarios, this effect is As shown in scenarios S8-S10, nitrate loads leaching from soil and transported by SGD are much more sensitive to land use changes than recharge and SGD flow rates. Our simulations show that converting forests and meadows to agriculture (S8) results in nitrate loads increased by about 75% with respect to the baseline scenario. This is consistent with the increase in agricultural area from 51% to 91% of the watershed. The large nitrate load is due to growing corn and canola, as part of crop rotation, assumed in the scenario.
Conversely, if crops and forests are replaced by grassland (S9), the nitrate loads are reduced by about six times compared to the baseline scenario. Grasslands are characterized by a high uptake of nitrogen, promoted by cutting hay twice a year, and the fertilized doses are also smaller.
The introduction of forests to the whole area (S10) also leads to the reduction of nitrate loads, albeit much less significant than in scenario S9. This result is caused by a smaller uptake of nitrogen in the grown-up forest as compared to growing grass. Even though there is no fertilization in the forest, a significant amount of nitrogen is released from peat soils. In other scenarios, this effect is offset with higher yields of crops and hay computed for peat areas. Further research is needed to quantify the actual rate of NO 3 leaching from peat soil in the considered area.
The values of recharge and N-NO 3 loads obtained from SWAT show significant spatial variability. Figure 7 presents the spatial distribution of values for a month with relatively high recharge and nitrate load (March 2007). The amount of recharge depends on the soil type and land use. Larger values are observed in sandy soils, mostly in the south-west part of the area. Lowest values occur predominantly in peat-covered river valleys. There is no simple relationship between recharge and nitrate load. For example, in the forests in the SW part of the area, nitrate leaching is limited, despite the moderate-to-high recharge. In contrast, in the southernmost part of the model domain, peat soil with a high content of organic matter releases large amounts of nitrate from mineralization of organic nitrogen, even though the recharge rate is small. offset with higher yields of crops and hay computed for peat areas. Further research is needed to quantify the actual rate of NO3 leaching from peat soil in the considered area. The values of recharge and N-NO3 loads obtained from SWAT show significant spatial variability. Figure 7 presents the spatial distribution of values for a month with relatively high recharge and nitrate load (March 2007). The amount of recharge depends on the soil type and land use. Larger values are observed in sandy soils, mostly in the south-west part of the area. Lowest values occur predominantly in peat-covered river valleys. There is no simple relationship between recharge and nitrate load. For example, in the forests in the SW part of the area, nitrate leaching is limited, despite the moderate-to-high recharge. In contrast, in the southernmost part of the model domain, peat soil with a high content of organic matter releases large amounts of nitrate from mineralization of organic nitrogen, even though the recharge rate is small.

Conclusions
SWAT, MODFLOW-NWT, and MT3DMS computer codes were sequentially coupled, in order to develop a model of transient groundwater flow and nitrate transport in a coastal multi-aquifer hydrosystem. The SWAT model allowed the capture of effects of hypothetical changes in land use and crop type on groundwater fluxes (including SGD) and the associated nitrate loads. The main findings can be summarized as follows:  Groundwater recharge, SGD, and the corresponding nitrate loads show a distinct time variable pattern, with maximum recharge rates and NO3 leaching in late winter/early spring.  The average values of recharge and SGD fluxes are influenced more significantly by crop type grown on farmlands than by the changes in land use. The maximum relative difference between the 10 y average of SGD flux between different scenarios did not exceed 12%. In contrast, nitrate leaching from soil and nitrate transport via SGD shows a larger variability, strongly depending on crop type and land use.  The lowest N-NO3 load in SGD occurred for the hypothetical scenario with all land converted to grassland, and it was three times smaller than the largest load, corresponding to converting all land to growing crops. While we believe that the general trends observed in our simulations are realistic, we are also aware of the limitations of the presented approach. Sequential coupling applied in this work does not allow for a consistent representation of aquifer-stream interactions in SWAT and MODFLOW models. In addition, the simplified representation of water and nitrate movement in the vadose zone between the root zone and groundwater table used by SWAT may lead to inaccuracies in the calculated time variability of contaminant transport, especially in complex geological conditions of the investigated area [69][70][71]. Further research will focus on an improved description of these processes in our model.

Conclusions
SWAT, MODFLOW-NWT, and MT3DMS computer codes were sequentially coupled, in order to develop a model of transient groundwater flow and nitrate transport in a coastal multi-aquifer hydrosystem. The SWAT model allowed the capture of effects of hypothetical changes in land use and crop type on groundwater fluxes (including SGD) and the associated nitrate loads. The main findings can be summarized as follows:

•
Groundwater recharge, SGD, and the corresponding nitrate loads show a distinct time variable pattern, with maximum recharge rates and NO 3 leaching in late winter/early spring.

•
The average values of recharge and SGD fluxes are influenced more significantly by crop type grown on farmlands than by the changes in land use. The maximum relative difference between the 10 y average of SGD flux between different scenarios did not exceed 12%. In contrast, nitrate leaching from soil and nitrate transport via SGD shows a larger variability, strongly depending on crop type and land use.

•
The lowest N-NO 3 load in SGD occurred for the hypothetical scenario with all land converted to grassland, and it was three times smaller than the largest load, corresponding to converting all land to growing crops.
While we believe that the general trends observed in our simulations are realistic, we are also aware of the limitations of the presented approach. Sequential coupling applied in this work does not allow for a consistent representation of aquifer-stream interactions in SWAT and MODFLOW models. In addition, the simplified representation of water and nitrate movement in the vadose zone between the root zone and groundwater table used by SWAT may lead to inaccuracies in the calculated time variability of contaminant transport, especially in complex geological conditions of the investigated area [69][70][71]. Further research will focus on an improved description of these processes in our model.

Conflicts of Interest:
The authors declare no conflict of interest.