Characterizing Groundwater Interaction with Lakes and Wetlands Using GIS Modeling and Natural Water Quality Measurements

Wetlands provide many benefits, including flood attenuation, groundwater recharge, water-quality improvement, and habitat for wildlife. As their structure and functions are sensitive to changes in hydrology, characterizing the water budgets of wetlands is crucial to effective management and conservation. The groundwater component of a budget, which often controls resiliency and water quality, is difficult to estimate and can be costly, time-consuming, and invasive. This study used a GIS approach using a digital elevation model (DEM) and the elevations of lakes, wetlands, streams, and hydric soils to produce a water-table surface raster for a portion of the Itasca Moraine, Minnesota, U.S. The water-table surface was used to delineate groundwatersheds and groundwater flow paths for lakes and wetlands, and map recharge and discharge rates across the landscape. Specific conductance and pH, which depend on the hydrological processes that dominate a wetlands water budget, were measured in the field to verify this modeling technique. While the pH of surface waters varied in the study area, specific conductance increased from 16.7 to 357.5 μS/cm downgradient along groundwater flow paths, suggesting increased groundwater interaction. Our results indicate that basic GIS tools and often freely available public-domain elevation datasets can be used to map and characterize the interaction of groundwater in the water budgets of lakes and wetlands, as exemplified by the Itasca Moraine region. Combining this with grid cell-by-cell water balance provides a means to estimate recharge and discharge, thereby affording a way to quantify groundwater contribution to and from lakes and wetlands. Applied elsewhere, this cost-efficient technique can be used to assess the vulnerability of lakes and wetlands to changes in land use, groundwater development, and climate change.


Introduction
Wetlands provide benefits to society, such as flood attenuation, groundwater recharge, water-quality improvement, and habitat for wildlife, including those that are threatened or endangered [1]. Unfortunately, the benefits of wetlands were not adequately understood and the legal protection of wetlands was not a priority in the United States until the 1970s. As a result, it is estimated that the United States has lost approximately 50% of its wetlands since European settlement [1].
Characterizing the water budgets of lakes and wetlands is crucial for their effective management. It is especially important to understand the influence of groundwater in these systems as surface water and groundwater are interconnected and wetlands are sensitive to changes in their hydrology [2,3]. Groundwater influence can vary substantially among even adjacent and nearby wetlands, and often groundwater catchment areas, referred to as groundwatersheds [4], do not coincide with surface watersheds, complicating matters [5]. Characterizing the groundwater component of the wetland water budget can help natural resource managers better understand how these features will respond to changes in their hydrology and help assure effective wetland conservation.
Traditional techniques of water budget estimation require extensive equipment, longterm collection of field data, and labor for characterization, which can quickly become costly. In addition, some field areas may be difficult to access, and field visits and the installation of equipment may damage the ecosystem. Although numerical groundwater flow models provide a framework for thoroughly understanding wetland water budgets, model calibration and application generally requires the same or even more extensive field data (e.g., [6]). In this report, we demonstrate the use of GIS tools and the measurement of pH and specific conductance as a novel approach and method to characterize the contribution of groundwater in the water budgets of wetlands in a more cost-effective and less invasive way. As an example, we applied the technique within northern Minnesota's Itasca Moraine, where the interaction and contribution of groundwater in the water budgets of lakes and wetlands was characterized with a spatial modeling technique using several GIS tools. Specific conductance and pH of surface waters, which depend on the hydrological processes that dominate a wetland's water budget, were measured in the field and used to verify and support our spatial-modeling technique.

Background and Previous Work
Surface water and groundwater are interconnected in most types of terrain; effective water-resource management requires an understanding of their interaction [2]. For lakes and wetlands in glacial terrain, groundwater interacts in one of three ways ( Figure 1): through recharge (flow from surface water), discharge (flow from groundwater to surface water), or a combination where recharge and discharge will occur in different locations within the same surface water body (flow-through). Although this interaction commonly corresponds with topographic position, it depends largely on where the lake or wetland lies within the larger-scale groundwater-flow system [3].
Water 2021, 13, x FOR PEER REVIEW 2 of 31 among even adjacent and nearby wetlands, and often groundwater catchment areas, referred to as groundwatersheds [4], do not coincide with surface watersheds, complicating matters [5]. Characterizing the groundwater component of the wetland water budget can help natural resource managers better understand how these features will respond to changes in their hydrology and help assure effective wetland conservation. Traditional techniques of water budget estimation require extensive equipment, long-term collection of field data, and labor for characterization, which can quickly become costly. In addition, some field areas may be difficult to access, and field visits and the installation of equipment may damage the ecosystem. Although numerical groundwater flow models provide a framework for thoroughly understanding wetland water budgets, model calibration and application generally requires the same or even more extensive field data (e.g., [6]). In this report, we demonstrate the use of GIS tools and the measurement of pH and specific conductance as a novel approach and method to characterize the contribution of groundwater in the water budgets of wetlands in a more costeffective and less invasive way. As an example, we applied the technique within northern Minnesota's Itasca Moraine, where the interaction and contribution of groundwater in the water budgets of lakes and wetlands was characterized with a spatial modeling technique using several GIS tools. Specific conductance and pH of surface waters, which depend on the hydrological processes that dominate a wetland's water budget, were measured in the field and used to verify and support our spatial-modeling technique.

Background and Previous Work
Surface water and groundwater are interconnected in most types of terrain; effective water-resource management requires an understanding of their interaction [2]. For lakes and wetlands in glacial terrain, groundwater interacts in one of three ways ( Figure 1): through recharge (flow from surface water), discharge (flow from groundwater to surface water), or a combination where recharge and discharge will occur in different locations within the same surface water body (flow-through). Although this interaction commonly corresponds with topographic position, it depends largely on where the lake or wetland lies within the larger-scale groundwater-flow system [3]. For areas with minimal subsurface information but abundant surface water features, a digital elevation model (DEM) can be used to roughly map a generalized distribution of groundwater recharge and discharge [7] based simply on topography. DEMs allow the delineation of topographic (surface) watersheds and provide the elevations of surface water bodies such as lakes, streams, and wetlands, which can serve as a basis for interpolating the surrounding water-table surface. Interpolation tools available in GIS programs can also be used to generate a water-table surface from the elevations of surface water features derived from a DEM.
The DEM resolution, or cell size, has been shown to have an effect on hydrological interpretations derived from topographic features [8][9][10]. Slope, upslope area, flow For areas with minimal subsurface information but abundant surface water features, a digital elevation model (DEM) can be used to roughly map a generalized distribution of groundwater recharge and discharge [7] based simply on topography. DEMs allow the delineation of topographic (surface) watersheds and provide the elevations of surface water bodies such as lakes, streams, and wetlands, which can serve as a basis for interpolating the surrounding water-table surface. Interpolation tools available in GIS programs can also be used to generate a water-table surface from the elevations of surface water features derived from a DEM.
The DEM resolution, or cell size, has been shown to have an effect on hydrological interpretations derived from topographic features [8][9][10]. Slope, upslope area, flow lengths, and the watershed area are commonly used in hydrological studies; these parameters are Water 2021, 13, 983 3 of 30 sensitive to DEM resolution [8,9,11]. For example, Wu et al. [9] found that mean slope decreased with decreasing DEM resolution and upslope area increased with decreasing DEM resolution. Zhang and Montgomery [8] found that coarser resolution DEMs led to a decrease in depth to the water table and an increase in peak discharge. These observations emphasize the importance of considering the sensitivity of a model to the DEM resolution and selecting an appropriate resolution for modeling objectives.
Water quality has also been found to be useful in groundwater modeling. Winter and Carr [12] created a groundwater model for moraine sediments in Missouri Coteau, North Dakota, and measured water quality to evaluate how lakes and wetlands interact with groundwater. Lakes with large quantities of total dissolved solids (TDS) received discharge from regional or deep flow systems while lakes with smaller amounts of TDS received discharge from smaller, more localized flow systems [2]. Within the Missouri Coteau study area, Winter and Carr [12] mapped recharge, discharge, and flow-through lakes and wetlands based on data from 1979.
More recently, Jaworksa-Szulc [13] used a groundwater model and TDS measurements to characterize how lakes interact with groundwater in the Young Glacial Area of northern Poland. Lakes that received groundwater discharge had TDS measurements that ranged from 250-350 mg/L while losing lakes or those that recharge groundwater measured less than 100 mg/L. The flow-through lakes identified had TDS measurements ranging from 170-200 mg/L. In addition, the discharge lakes were commonly found in ancient tunnel valleys or other meltwater-derived features while the losing (recharge) lakes were found on ground moraine.
Gerla [14] suggested that pH and electrical conductivity (or specific conductance), both easily measured in the field, can be useful in determining what processes are dominant in a wetland's water budget. Groundwater interacts with subsurface materials and dissolves minerals as it flows, thereby increasing its specific conductance. Generally, groundwater will have a much larger specific conductance compared to surface water [15]. If the water budget is dominated by meteoric water (ombrotrophic), then the waters of the wetland will be less mineralized and have a characteristically low specific conductance. If the water budget is dominated by groundwater (minerotrophic), the wetland waters will be more mineralized and have an elevated specific conductance [2,[16][17][18].
Meteoric water interacts with carbon dioxide, other gases, and particulates in the atmosphere, which results in a slightly acidic pH of approximately 5.6 [19,20]. Due to this general trend, wetlands with water budgets dominated by meteoric water, such as ombrotrophic bogs, are more acidic [17,18]. Several factors affect surface water and groundwater pH. For example, aquatic plant productivity in surface waters influences pH as photosynthetic reactions consume carbon dioxide, resulting in an alkaline pH [21,22]. In bogs, vegetation such as sphagnum moss has the opposite effect and maintains acidity [23][24][25]. This process is enhanced by limited groundwater seepage which would otherwise introduce neutralizing bases [21]. Water in wetlands that is contributed from groundwater is also susceptible to pH variation, depending on the composition of the soil, sediment, or rock with which the water has interacted (e.g., [26]). It also is influenced by subsurface residence time, the contributing basin size, and biological activity [3,27,28].
These observational studies of TDS and pH show that basic water quality data may be useful in delineating groundwater-surface water interaction and flow paths, although these prior studies do not couple this tool together with water balance, topographic analysis, and/or spatial modeling.
Our work adds to existing knowledge by providing a simple, novel spatial modeling approach that can delineate surface water-groundwater interaction. This technique can be applied in areas with abundant surface water features and is low-cost and less invasive than traditional techniques used to characterize surface water and groundwater interaction. This research focuses on the interaction between groundwater and lakes and wetlands in the Itasca Moraine and provides insight into the value of using pH and specific conductance as a tool for verifying our technique.

Study Site
We selected the Elk Lake and Nicolet Creek watersheds for spatial modeling. These contiguous watersheds together cover 47 km 2 and fall mostly within Minnesota's Itasca State Park, named after Lake Itasca and proclaimed by Henry Schoolcraft in 1832 as the headwaters of the Mississippi River [29]. The Nicolet Creek and Elk Lake watersheds, which lie entirely within the Itasca Moraine, drain directly into Lake Itasca, and host many small lakes, wetlands, creeks, and groundwater springs. Unlike many other areas in the region, these watersheds have been protected from agricultural drainage and changes to the original forest land cover. Due to the landscape creation during the Late Pleistocene, surface drainage in the watersheds is poorly integrated and characterized by shallow groundwater flow. The water table is exposed in lakes, ponds, and wetlands that cover roughly 25-30% of the area, which create the primary hydrological data for our model. The watersheds lie at elevations that range from 447 to 550 m-the highest in the region, making it highly unlikely that any groundwater from any distant area is discharged in either watershed. Few hydrological data are available for the watersheds, except for Elk Lake, which has more than 80 years of water level records. During that period (1938-2020), the maximum variation was only 0.94 m (https://www.dnr.state.mn.us/lakefind/lake.html?id=15001000, accessed on 27 December 2020), indicating a strongly static water level regime. The region's midcontinental climate is characterized by warm humid summers and cold winters, with most of the 0.75 m yearly average precipitation occurring during convective summer storms.
The Itasca Moraine ( Figure 2) was deposited by the stagnation of the Wadena Lobe and was subsequently altered by regional glacial activity until the end of the Wisconsinan glaciation, approximately 11,000 years ago [30,31]. The Itasca Moraine exhibits approximately 200 m of topographic relief and extends for nearly 150 km across northern Minnesota where it serves as a regional recharge zone ( Figure 2). The Wadena lobe till that composes the moraine is rich in Paleozoic carbonate fragments derived glacially from southern Manitoba, Canada. It also contains Precambrian igneous and metamorphic fragments which are believed to have been eroded from the Canadian Shield. The water chemistry of Elk Lake is largely influenced by the calcareous glacial drift of the moraine, with Ca 2+ , Mg 2+ , and HCO 3 − as dominant ions. The lake precipitates CaCO 3 and its varved sediments have been studied extensively to elucidate changes in Holocene climate [32].
The study area is underlain by end moraine and outwash deposits [33] (Figure 2) and has a hummocky topography characterized by numerous small lakes and wetlands that were formed when ice blocks broke off from the ice lobe and later melted, leaving kettle depressions [30]. Larger lakes, such as Lake Itasca and Elk Lake, occupy former subglacial channels, locally referred to as tunnel valleys, which were likely formed by drainage within the ablating ice sheet [30]. Landform assemblages for the study area were described by [31] and consist of outwash channels and a supraglacial complex. The outwash channels are primarily well-sorted sand and gravel and the supraglacial complex is primarily sandy loam that is poorly to moderately sorted.  The map on the right shows the Elk Lake and Nicolet Creek (topographic/surface watershed outlined in black). Surficial Quaternary geology map of the study area is from Hobbs and Goebel (1982) overlain with lake and wetland areas (light blue), which were used to interpolate the water-table surface.
The study area is underlain by end moraine and outwash deposits [33] (Figure 2) and has a hummocky topography characterized by numerous small lakes and wetlands that were formed when ice blocks broke off from the ice lobe and later melted, leaving kettle depressions [30]. Larger lakes, such as Lake Itasca and Elk Lake, occupy former subglacial channels, locally referred to as tunnel valleys, which were likely formed by drainage within the ablating ice sheet [30]. Landform assemblages for the study area were described by [31] and consist of outwash channels and a supraglacial complex. The outwash channels are primarily well-sorted sand and gravel and the supraglacial complex is primarily sandy loam that is poorly to moderately sorted.

Approach and Data
The work flow of our method is shown in Figure 3 and described in detail in the following subsections. There are four overall steps in the GIS-based technique we used for investigating groundwater-surface water interaction: (1) Generate a local and regional water-table surface from a DEM and associated lake, stream, wetland, and hydric soil elevation data. (2) Use the water-table surface to delineate groundwatersheds for lakes and wetlands (blue arrows), which are then compared to topographic (surface) watersheds (green arrows) to determine the relative importance of groundwater and surface water for selected lakes.  The map on the right shows the Elk Lake and Nicolet Creek (topographic/surface watershed outlined in black). Surficial Quaternary geology map of the study area is from Hobbs and Goebel (1982) overlain with lake and wetland areas (light blue), which were used to interpolate the water-table surface.

Approach and Data
The work flow of our method is shown in Figure 3 and described in detail in the following subsections. There are four overall steps in the GIS-based technique we used for investigating groundwater-surface water interaction: (1) Generate a local and regional water-table surface from a DEM and associated lake, stream, wetland, and hydric soil elevation data. (2) Use the water-table surface to delineate groundwatersheds for lakes and wetlands (blue arrows), which are then compared to topographic (surface) watersheds (green arrows) to determine the relative importance of groundwater and surface water for selected lakes. (3) Based on the water-table surface, flow accumulation is computed. Groundwater flow paths are derived from this step and coupled with the location of lakes and wetlands. (4) Recharge-discharge across the landscape is mapped by combining data from the water-table surface, saturated soil hydraulic conductivity, and an analytical technique to estimate the maximum depth of groundwater-surface water interaction.
Lastly, the output shows the strength and proximity of recharge and discharge areas relative to lakes and wetlands (red arrows), which are then corroborated by specific conductance and pH measurements collected in the field. Lastly, the output shows the strength and proximity of recharge and discharge areas relative to lakes and wetlands (red arrows), which are then corroborated by specific conductance and pH measurements collected in the field. Public domain datasets were used for the groundwater model, including a 3-m DEM derived from LiDAR (light ranging and detection) and obtained from the Minnesota Geospatial Information Office's MnTOPO Viewer [34] ( Figure S1). Lake and wetland data were obtained from the Minnesota National Wetland Inventory Update [35], which is based on high-resolution aerial imagery and LiDAR, creating the most comprehensive, current, and accurate wetland inventory in the U.S. The northwest region, which includes the Elk Lake and Nicolet Creek watersheds, was released to the public in January 2019. Summer and fall 2017 statewide aerial imagery (1 m resolution), provided by the Minnesota Geospatial Information Office's Web Mapping Service [36], was used to identify and verify lakes and wetlands. A shapefile including streamlines was obtained from the National Hydrography Dataset (NHD) [37]. Soil polygons (1:20,000 and 1:24,000 scale) for the study area came from the Soil Survey Geographic Database (SSURGO) [38]. The SSURGO soil data table was used to estimate hydraulic conductivity and produce the recharge-discharge map [38,39].

Interpolation of the Water-Table Surface
Lakes were the primary source of elevation data for deriving the water-table surface. Wetlands with specific hydrologic modifier codes that indicate permanent saturation were included to supplement the lake data. Hydrologic modifier codes of wetlands selected for inclusion were those classified as F-semipermanently flooded, G-intermittently exposed, and H-permanently flooded using the Cowardin et al. system [40]. All other wetlands were removed from the lake and wetland dataset. An additional raster dataset that mapped the minimum elevation value for each wetland and lake polygon was created to eliminate the small range of elevation variability within the DEM. This ensured that the lowest topographic elevation present within each polygon was assigned to the entire water body. Public domain datasets were used for the groundwater model, including a 3-m DEM derived from LiDAR (light ranging and detection) and obtained from the Minnesota Geospatial Information Office's MnTOPO Viewer [34] ( Figure S1). Lake and wetland data were obtained from the Minnesota National Wetland Inventory Update [35], which is based on high-resolution aerial imagery and LiDAR, creating the most comprehensive, current, and accurate wetland inventory in the U.S. The northwest region, which includes the Elk Lake and Nicolet Creek watersheds, was released to the public in January 2019. Summer and fall 2017 statewide aerial imagery (1 m resolution), provided by the Minnesota Geospatial Information Office's Web Mapping Service [36], was used to identify and verify lakes and wetlands. A shapefile including streamlines was obtained from the National Hydrography Dataset (NHD) [37]. Soil polygons (1:20,000 and 1:24,000 scale) for the study area came from the Soil Survey Geographic Database (SSURGO) [38]. The SSURGO soil data table was used to estimate hydraulic conductivity and produce the recharge-discharge map [38,39].

Interpolation of the Water-Table Surface
Lakes were the primary source of elevation data for deriving the water-table surface. Wetlands with specific hydrologic modifier codes that indicate permanent saturation were included to supplement the lake data. Hydrologic modifier codes of wetlands selected for inclusion were those classified as F-semipermanently flooded, G-intermittently exposed, and H-permanently flooded using the Cowardin et al. system [40]. All other wetlands were removed from the lake and wetland dataset. An additional raster dataset that mapped the minimum elevation value for each wetland and lake polygon was created to eliminate the small range of elevation variability within the DEM. This ensured that the lowest topographic elevation present within each polygon was assigned to the entire water body.
Streams data obtained from the NHD [37] were used to add elevation control and for gridding surface breaklines. Stream lines were inspected and edited to ensure that they did not transect lakes and wetlands. Most stream lines were edited to include more vertices or redrawn to more closely match stream channels observed in the aerial imagery.
SSURGO hydric soil maps [38] were also used to gain greater elevation control and add constraint on the modeled water-table elevation. A hydric soil is one that has developed anaerobic conditions and features as a result of saturation during the growing season [41]. Soil polygons were included if they met the following criteria: 1.
The polygon was listed as a hydric soil and was described as frequently ponded, which is a closed depression that contains water and has more than a 50% chance of occurring in any given year under normal weather conditions [41]; 2.
The soils represented by the polygon covered small parts of the raster where additional water-table elevation control was necessary, including areas where the water-table surface lay above the DEM.
Iterative finite difference interpolation [42] was used to create the water-table raster from the lake, wetland, and stream data using GIS tools. For comparison, this provisional water-table surface raster was subtracted from the DEM (Figure 3), the difference showing areas where the water table lies above the DEM, indicating the input needed adjustment ( Figure S2). In the few scattered areas where this occurred, polygons were constructed in areas where standing water was observed and stream lines created where dendritic drainage features could be discerned on aerial imagery. These features were added and the entire dataset reinterpolated (diamond symbol, Figure 3).
Interpolating the water-table surface was an iterative process. Only two to five additional features were added for elevation control at a time before re-running the interpolation tool. The resulting water-table surface was then again compared to the DEM to determine what additional improvements were made to the resulting surface. The process took approximately 30 trials using the interpolation tool to create the final water-table surface. In total, 883 polygons were used to interpolate the water-table surface, including 761 wetland inventory polygons, 102 hydric soil polygons, and 20 polygons created from the aerial imagery and the DEM ( Figure S3).
Interpolation was carried out using ANUDEM, an algorithm originally developed by the Australia National University [42] for the construction of hydrologically correct DEMs (e.g., [43]). This algorithm is essentially a bivariate thin plate smoothing spline that uses finite difference discretization. It is more efficient and maintains a high-quality surface continuity when compared to other global interpolation methods, including thin plate spline and kriging when using large datasets [42,44].
When this technique is applied to create a water-table surface, the algorithm uses the known elevations of the input surface water features to estimate the elevation of the water table between these features. The algorithm assigns elevation data supplied by the surface water features to the nearest grid point. Gauss-Seidel iteration with overrelaxation [45,46] is used to calculate nearby unknown grid points. The process starts with a coarse grid and the spacing is diminished sequentially until the user-specified grid resolution is obtained. A drainage enforcement algorithm operates simultaneously and works to remove sinks in the surface through imposing drainage conditions [42,44].
To create the water-table surface, the lake and wetland elevation data were input as both contours and lakes. The input contours are used at first to determine the general morphology of the water-table surface using points of maximum curvature. The contours are then used to estimate the elevations of cells in between input features and the generalized morphology is updated iteratively throughout this process [42,44]. Inputting the lake and wetland polygons as lakes ensured that the elevation of the interior cells lie below the polygon boundary while the elevation of the outside cells lie above it [42,44]. The input of stream data assists the interpolation process by serving as a gridding break line and ensures the removal of data that conflict with the streamlines' descent [42,44]. The elevations of the stream line vertices are also assigned to the nearest grid point to assist with interpolation [44].

Groundwater Flow Pathlines and Groundwatersheds
The water-table surface was used to generate groundwater-flow pathlines using techniques that are well established for topographic DEM channels [47]. To produce a connected drainage network, it was necessary to fill sinks present in the water-table surface. Sinks often result from the rounding of elevation values when generating a surface but can also represent actual depressions in the landscape. In either case, these features result in a discontinuous drainage network if they are not filled [47,48].
An eight-direction (D8) flow model was then used to determine the direction of flow out of each cell, which is determined by the steepest descent [48]. This information was then used to determine the flow accumulation, or the number of cells that flow into a particular cell. Cells with high accumulation are areas were flow is concentrated and can be used as flow paths [47,48]. Cells that received accumulated flow from 250 or more cells defined the flow paths, which ensured that all the wetlands sampled had a flow path in their vicinity. The cumulative flow path length for each sampled water body was then calculated by summing the lengths of associated flow paths. If multiple flow paths converged at a lake or wetland, the length for the longest reach was summed. In the case of Elk Lake, which had two large tributaries, a weighted average was calculated using the contributing groundwater drainage area for each tributary.
The flow path lengths were then plotted against the average specific conductance and pH measurements for each water body for statistical comparison. Linear regression was used to analyze the strength of the relationship between groundwater flow path lengths and pH and specific conductance measured in the field. The strength of the relationship was evaluated with R 2 and the p-value. R 2 provided the percent of the variance of pH and specific conductance that is explained by the effect of the groundwater flow path length. The p-value provided the likelihood of obtaining water quality measurements when the null hypothesis is true [49]. In this case, the null hypothesis was that there is no relationship between the groundwater flow path length and either pH or specific conductance for lakes and wetlands. If the p-value was computed to be less than the selected significance value of 0.05, the results were assumed to be statistically significant and the null hypothesis was rejected.
Topographic (surface) watersheds and groundwatersheds were delineated using the same method, but with different datasets. Surface watersheds were based on the 3-m DEM directly, while the interpolated water-table surface raster was used for groundwatersheds. Both were delineated by identifying all upslope raster cells (eight-point flow direction (D8) model [48]) from an outlet point. The outlet for the groundwatershed was identified by locating the largest flow accumulation value for the lake or wetland of interest. In the case of a lake or wetland with no outlet stream, then the lake or wetland feature itself was selected as the outlet. Either case, a raster containing upslope cells defined the groundwatershed. Finally, for each water body sampled, its surface watershed to groundwatershed area ratio was computed and compared statistically to the average specific conductance and pH measurements. Linear regression was used to analyze the strength of the relationship between the log-transformed ratio computed from the surface watershed and groundwatershed area and the pH and specific conductance measured in the field. The strength of the relationship was evaluated with R 2 and considered significant if p < 0.05.

Mapping Recharge-Discharge
Darcy's Law describes flow through a porous medium and can be combined with the conservation of mass to determine water balance for cells within a raster. Equation (1) shows Darcy's Law, where groundwater discharge Q (m 3 /s) is a function of hydraulic conductivity K (m/s), the cross-sectional area perpendicular to flow A (m 2 ), and the hydraulic gradient ∆h L , where ∆h (m) is the difference in hydraulic head between cells and L (m) is the distance between cell centers (m): Conductance, C, is defined by Equation (2) [50]: where h A and h B are the hydraulic heads at the ends of a prism ( Figure 4a) that makes up one cell of the discretized groundwater flow system. Considering a stack of subprisms or group of cells along the direction of groundwater flow (Figure 4b), then Equation (3) gives the continuity of head across the cells: where i = 1 ton enumerates the subprism group of cells ( Figure 4b). Assuming one-dimensional flow and no change in storage, combining Equations (2) and (3) gives Equation (4): and applying Equation (2) with rearrangement gives: where c i represents the conductance of the ith subprism. Equation (5) implies that for a set of conductances arranged in a series, the inverse of the equivalent conductance equals the sum of the individual conductances, which is analogous with a harmonic mean. If there are only two cells (subprisms), then Equation (6) applies: where c 1 and c 2 are the conductances of the two cells. Thus, for two dimensional flow in a equilinear grid, Darcy's Law for heterogeneous confined flow in the x-axis direction can be determined by Equation (7): where b is the vertical thickness of the layer, i + 1/2 represents the flow across the i and i + 1 cell-to-cell boundary in the x-axis direction, and j enumerates the grid cells in the y-axis direction. A similar relationship can be written for flow in the y-axis direction. Groundwater flow in the Itasca Moraine is generally unconfined because of the thick, relatively coarse surficial glacial deposits [30] and the numerous lakes and wetlands assumed to expose the water table. For unconfined conditions, saturated thickness and hydraulic head are dependent. The Dupuit approximation [19,51] linearizes this dependency. Due to the gentle gradient of the water table and the lack of steep topographic slopes, the Dupuit approximation is likely valid for most of the Elk Lake and Nicolet Creek watersheds. Thus, the calculation of unconfined specific discharge q , or flow per unit width (length 2 /time) from hydraulic conductivity and the hydraulic gradient applies: The unconfined cell-to-cell discharge can then be approximately determined by combining Equations (7) and (8): where K i,j represent the hydraulic conductivities in the two adjacent cells along the x-axis direction. A similar equation can be written for flow in the y-axis direction, with cells enumerated by j. In two dimensions, Darcy's Law can be applied to each of the four sides of each raster cell and the four discharge values summed to give a flow residual R (m 3 /s) (Equation (10)) for the cell:  Groundwater flow in the Itasca Moraine is generally unconfined because of the thick, relatively coarse surficial glacial deposits [30] and the numerous lakes and wetlands assumed to expose the water table. For unconfined conditions, saturated thickness and hydraulic head are dependent. The Dupuit approximation [19,51] linearizes this dependency. Due to the gentle gradient of the water table and the lack of steep topographic slopes, the Dupuit approximation is likely valid for most of the Elk Lake and Nicolet Creek (c) Example of groundwater recharge, discharge, and balanced cells (shaded) in a two-dimensional volume balance residual raster, where conductance is homogeneous, h(i + 1,j) and h(i,j + 1) > h(i,j), and h(i − 1,j) and h(i,j − 1) < h(i,j). The bolder arrows represent greater flow magnitude.
Application of this procedure on the water-table map generated the groundwater balance residual raster across the flow domain in which the residuals are the differences between inflow and outflow of the cells (Figure 4c) and was used to create a recharge-discharge map. The imbalance of groundwater flow within the cell must be accommodated by either inflow (recharge) or outflow (discharge) to comply with the conservation of mass. A negative residual indicates that recharge from surface water exceeds discharge out of the cell while a positive residual indicates that discharge from the cell exceeds recharge from surface water.
Hydraulic conductivity was estimated from data reported in SSURGO [38] ( Figure S4). Soils in the database include both field and laboratory investigation of the soil's characteristics [39]. Wieczorek [39] used soil information from SSURGO, including Ksat, and produced an area and depth weighted table of soil characteristics that can be appended to soil polygons for any area of interest in the United States. This was done specifically for studies that require spatial estimates of soil hydraulic properties [39]. No Ksat values were available for lake beds within the study area, so they were set to be similar to that of wetlands in the study area, 15 µm/s or 1.3 m/day.
The interpolated water-table surface and an estimated aquifer base elevation were needed to create the hydraulic-head raster. In three transects of the Itasca Moraine within the park, we applied an analytical model developed by Bresciani et al. [52], which is based on Hooghoudt's drainage equation [52,53]. Data required include saturated hydraulic conductivity, based on SSURGO [38], recharge estimates [54], and water-level elevations from the 3-m DEM. The analytical model (Appendix A) provided a best for fit surface water and groundwater interaction at various depths but reached an optimal value at 441 m (above mean sea level) for the transects that crossed the study area. This elevation was assumed to be the base of the unconfined aquifer in subsequent models. The unconfined hydraulic-head raster was created by subtracting the aquifer base elevation, 441 m, from the interpolated water-table surface.
As hydrological interpretations have been shown to be sensitive to the resolution of elevation rasters [8,9,11], the sensitivity of the recharge-discharge map to the water-table surface cell size was determined. This was accomplished by summing the water balance for all cells within the groundwatershed (Appendix A) and then comparing it to streamflow measurements at Nicolet Creek and the Elk Lake outlet, Chambers Creek. This approach ensures that the water balance for the groundwater volume balance residual raster is similar to the discharge leaving the two outlets. In an effort to better match the discharge measurements of the two outlets, the cell size for the water-table surface raster was varied. Sizes selected for the sensitivity analysis were 10, 20, 30, 40, and 50 m.

Field Methods
The wetland dataset discussed previously was used to select wetlands for sampling. A Hanna HI9813-6 pH/SC/TDS/temperature meter (Ann Arbor, MI, USA) was used to measure pH and specific conductance in approximately 20 groundwatershed wetlands in September 2018. The meter was calibrated prior to the start of fieldwork using a 7.01 pH buffer and a 1413 µS/cm conductivity standard. Ultimately, lakes and wetlands sampled were selected by their Cowardin et al. [40] hydrologic modifier code (F-semi-permanently flooded, G-intermittently exposed, and H-permanently flooded), location within the watershed, and accessibility. Effort was made to select an adequate distribution of lakes and wetlands across the study area and a small kayak was used so that measurements could be taken in a variety of locations within each waterbody. Measurements were taken at 5-10 locations in each waterbody, depending on its size. This included at least one measurement in the center and one in each quadrant ( Figure 5).
buffer and a 1413 μS/cm conductivity standard. Ultimately, lakes and wetlands sampled were selected by their Cowardin et al. [40] hydrologic modifier code (F-semi-permanently flooded, G-intermittently exposed, and H-permanently flooded), location within the watershed, and accessibility. Effort was made to select an adequate distribution of lakes and wetlands across the study area and a small kayak was used so that measurements could be taken in a variety of locations within each waterbody. Measurements were taken at 5-10 locations in each waterbody, depending on its size. This included at least one measurement in the center and one in each quadrant ( Figure 5). Figure 5. Wetland, pond, spring, and stream field measurement locations discussed in the text. Features referred to in the text: 1, Unnamed NW Pond 1; 2, Unnamed NW Pond 2; 3, Unnamed NW Pond 3; 4, Bog near Whipple Lake; 5, Nicolet Lake; 6, Triplet Lakes; 7, Unnamed Pond near Little Elk Lake; 8, Unnamed Pond near Whipple Lake; 9, Whipple Lake; 10, Augusta Lake; 11, Unnamed SW Pond 1; 12, Unnamed SW Pond 2; 13, Unnamed SW Pond 3; 14, Horn Lake; 15, Hernando DeSoto Lake; 16, Nicolet Creek; 17, Clark Lake; 18, Unnamed Pond near Clark Lake; 19, Nicolet (Chambers) Creek; 20, Elk Lake Spring 1; 21, Elk Lake Spring; 22, Elk Lake Spring 3; 23, Elk Lake Spring 4; and 24, Elk Lake.
The following steps were taken when conducting field measurements using published Minnesota Pollution Control Agency protocol [55] and the meter operation manual [56]: 1. The probe was rinsed with deionized water before and after each measurement; 2. Algae, floating vegetation, and disturbance of bottom sediments were avoided when taking measurements; 3. For each measurement of pH and specific conductance, the probe was placed approximately 10 cm below the water surface and swirled gently until the values stabilized, approximately 60 s; Figure 5. Wetland, pond, spring, and stream field measurement locations discussed in the text. The following steps were taken when conducting field measurements using published Minnesota Pollution Control Agency protocol [55] and the meter operation manual [56]: 1.
The probe was rinsed with deionized water before and after each measurement; 2.
Algae, floating vegetation, and disturbance of bottom sediments were avoided when taking measurements; 3.
For each measurement of pH and specific conductance, the probe was placed approximately 10 cm below the water surface and swirled gently until the values stabilized, approximately 60 s; 4.
At least one measurement was taken near the center of the waterbody and additional measurements were taken in each quadrant, at each site coordinates were recorded for each location with a Garmin eTrex 30x GPS receiver (Olathe, KS, USA).
Measurements at Elk Lake were taken from along the shore and at a dock. Streams and springs in the study area were also measured for comparison using methods similar to that of the lakes. For springs, the probe was submerged into the flowing spring water, without disturbing underlying sediments.

Water-Table Surface
In the study area, the water-table elevation ranges from approximately 448 m near Elk Lake to its maximum 520 m at the far west end of the study area. The groundwatershed for Nicolet Creek at its outlet on Lake Itasca and Elk Lake was delineated using the interpolated water-table surface ( Figure 6).
The best-fit interpolated water-table surface was subtracted from the DEM to map the depth to the water table, which shows a few scattered areas where the interpolated water-table surface lies above the surface topography ( Figure S5), although it occurred in less than 2% of the groundwatershed raster cells. Of these cells, 97.5% were less than 3 m above the DEM. No cells exceeded 5 m above the DEM.

Water-Table Surface
In the study area, the water-table elevation ranges from approximately 448 m near Elk Lake to its maximum 520 m at the far west end of the study area. The groundwatershed for Nicolet Creek at its outlet on Lake Itasca and Elk Lake was delineated using the interpolated water-table surface ( Figure 6). The best-fit interpolated water-table surface was subtracted from the DEM to map the depth to the water table, which shows a few scattered areas where the interpolated water-table surface lies above the surface topography ( Figure S5), although it occurred in less than 2% of the groundwatershed raster cells. Of these cells, 97.5% were less than 3 m above the DEM. No cells exceeded 5 m above the DEM.
The hydraulic gradient shows several areas with values greater than 0.05 on the north end of the study area (Figure 7). These elongate, north-south oriented areas range up to 3 km long, with the largest feature lying to the west of Nicolet Creek. Several smaller features lie to the east and southeast of Elk Lake with maximum gradients that also exceed 0.05. These features appear to correspond to mapped contacts between end moraine and glacial outwash sediments that extend north-south through the study area [33]. Locations of groundwater springs from the Minnesota Spring Inventory [57] and springs that were sampled in the study area generally coincide with areas at the base of large water-table slopes ( Figure 7). The hydraulic gradient shows several areas with values greater than 0.05 on the north end of the study area ( Figure 7). These elongate, north-south oriented areas range up to 3 km long, with the largest feature lying to the west of Nicolet Creek. Several smaller features lie to the east and southeast of Elk Lake with maximum gradients that also exceed 0.05. These features appear to correspond to mapped contacts between end moraine and glacial outwash sediments that extend north-south through the study area [33]. Locations of groundwater springs from the Minnesota Spring Inventory [57] and springs that were sampled in the study area generally coincide with areas at the base of large water-table slopes (Figure 7).

Groundwater Flow Paths
Calculated path lines of groundwater flow (Table 1) show northward flow with convergence towards Elk Lake and Lake Itasca (Figure 8). Ranging from 6.7 to 8.9, pH varied greatly but did not show any statistically significant relationship to the length of the up-

Groundwater Flow Paths
Calculated path lines of groundwater flow (Table 1) show northward flow with convergence towards Elk Lake and Lake Itasca (Figure 8). Ranging from 6.7 to 8.9, pH varied greatly but did not show any statistically significant relationship to the length of the upgradient groundwater flow path (R 2 = 0.18, p = 0.12) ( Figure S6).  Lengths of groundwater flow paths were also compared with average specific conductance measurements (Figure 9) and showed a statistically significant although weak relationship (R 2 = 0.39, p = 0.01) (Figure 10). Specific conductance appears to be larger in lakes and wetlands that have longer upgradient flow paths, although there are a few ex- Lengths of groundwater flow paths were also compared with average specific conductance measurements (Figure 9) and showed a statistically significant although weak relationship (R 2 = 0.39, p = 0.01) ( Figure 10). Specific conductance appears to be larger in lakes and wetlands that have longer upgradient flow paths, although there are a few exceptions. A few lakes (e.g., Augusta Lake and Horn Lake) that lie near the beginning of flow path lines have a relatively large specific conductance. In addition, unnamed ponds in the northwest portion of the study area lie at the middle of flow paths but have the lowest specific conductance measurements sampled in the study area (Figure 9).

Watersheds
The surface watershed-to-groundwatershed area ratio (SW:GW) was estimated for the lakes and wetlands that were sampled for pH and specific conductance (Table 1). For example, the model results show that Elk Lake has a similarly sized surface watershed and groundwatershed, and a SW:GW ratio of 1.12 was estimated ( Figure 11). In contrast, the model showed that some of the ponds have a larger surface watershed compared to its groundwatershed (Figure 12a) and others have a groundwatershed that is estimated to be as much as 25 times larger than its surface watershed (Figure 12b). For the lakes and wetlands that were sampled, SW:GW ranged from 0.04 to 6.11 ( Figure S7).

Watersheds
The surface watershed-to-groundwatershed area ratio (SW:GW) was estimated for the lakes and wetlands that were sampled for pH and specific conductance (Table 1). For example, the model results show that Elk Lake has a similarly sized surface watershed and groundwatershed, and a SW:GW ratio of 1.12 was estimated ( Figure 11). In contrast, the model showed that some of the ponds have a larger surface watershed compared to its groundwatershed (Figure 12a) and others have a groundwatershed that is estimated to be as much as 25 times larger than its surface watershed (Figure 12b). For the lakes and wetlands that were sampled, SW:GW ranged from 0.04 to 6.11 ( Figure S7).

Watersheds
The surface watershed-to-groundwatershed area ratio (SW:GW) was estimated for the lakes and wetlands that were sampled for pH and specific conductance (Table 1). For example, the model results show that Elk Lake has a similarly sized surface watershed and groundwatershed, and a SW:GW ratio of 1.12 was estimated ( Figure 11). In contrast, the model showed that some of the ponds have a larger surface watershed compared to its groundwatershed ( Figure 12a) and others have a groundwatershed that is estimated to be as much as 25 times larger than its surface watershed (Figure 12b). For the lakes and wetlands that were sampled, SW:GW ranged from 0.04 to 6.11 ( Figure S7). Water 2021, 13, x FOR PEER REVIEW 17 of 31 Figure 11. Surface watershed (light blue outline) and groundwatershed (red outline) delineated for Elk Lake.   Linear regression was used to compare the log-transformed SW:GW to the pH and specific conductance field measurements. In both cases, results did not correlate (R 2 = 0.22 and R 2 = 0.13) (Figures S7 and S8). Note, however, the trend that suggests that specific conductance converges and becomes less variable as SW:GW increases. There also are a few small lakes with high pH measurements that skew the results. For example, if the Triplet Lakes pH datum is removed, the pH-log watershed ratio correlation becomes significant (R 2 = 0.50, p = 0.0022) indicating increasing pH with decreasing SW:GW.

Spatial Distribution of Recharge and Discharge
Recharge-discharge results show the residual flow after balancing cell-by-cell groundwater input and output ( Figure 13). Strong discharge occurs on the east side of Elk Lake and a prominent discharge area to the west of it and Nicolet Creek. These features appear to correspond with springs from the Minnesota Spring Inventory [57] and those sampled, as well as areas that were mapped as having a water-table gradient that exceeds 0.05. The springs measured in the study area exhibited the largest specific conductance measurements, averaging 580 µS/cm. In addition, one wetland sampled east of Whipple Lake mapped as a strong recharge area surrounded by discharge. The relatively low pH measured in the field (6.7) and the type of vegetation observed (Sphagnum sp., Picea mariana, Larix laricina, Betula papyrifera, among others) suggests that this wetland may be an ombrotrophic bog. Linear regression was used to compare the log-transformed SW:GW to the pH and specific conductance field measurements. In both cases, results did not correlate (R 2 = 0.22 and R 2 = 0.13) (Figures S7 and S8). Note, however, the trend that suggests that specific conductance converges and becomes less variable as SW:GW increases. There also are a few small lakes with high pH measurements that skew the results. For example, if the Triplet Lakes pH datum is removed, the pH-log watershed ratio correlation becomes significant (R 2 = 0.50, p = 0.0022) indicating increasing pH with decreasing SW:GW.

Spatial Distribution of Recharge and Discharge
Recharge-discharge results show the residual flow after balancing cell-by-cell groundwater input and output ( Figure 13). Strong discharge occurs on the east side of Elk Lake and a prominent discharge area to the west of it and Nicolet Creek. These features appear to correspond with springs from the Minnesota Spring Inventory [57] and those sampled, as well as areas that were mapped as having a water-table gradient that exceeds 0.05. The springs measured in the study area exhibited the largest specific conductance measurements, averaging 580 μS/cm. In addition, one wetland sampled east of Whipple Lake mapped as a strong recharge area surrounded by discharge. The relatively low pH measured in the field (6.7) and the type of vegetation observed (Sphagnum sp., Picea mariana, Larix laricina, Betula papyrifera, among others) suggests that this wetland may be an ombrotrophic bog. Figure 13. Recharge-discharge map for the study area computed from the interpolated water-table surface based on the unfilled 40 × 40-m resolution raster. Areas of bolder color indicate greater than regional recharge (orange) and discharge (magenta) estimates, lighter colors indicate less than average rates (based on regional analysis [54]).
As the raster scale affects the calculated magnitude of local groundwater flow paths, the cell size of the water-table surface raster influences the overall water balance for the Figure 13. Recharge-discharge map for the study area computed from the interpolated water-table surface based on the unfilled 40 × 40-m resolution raster. Areas of bolder color indicate greater than regional recharge (orange) and discharge (magenta) estimates, lighter colors indicate less than average rates (based on regional analysis [54]).
As the raster scale affects the calculated magnitude of local groundwater flow paths, the cell size of the water-table surface raster influences the overall water balance for the groundwatershed. When the sum of the raster cell residuals across the entire groundwatershed was compared to the surface water discharge measurements made at the outlets in the field (Table 2), a raster with a 40-m cell-size provided the best result. This consisted of 0.54 m 3 /s discharging from the model, compared to an average 0.38 m 3 /s discharge rate measured and summed for the two outlets (Tables 2 and 3).  Table Surface Abundant lakes, wetlands, and streams, along with the availability of a 3-m DEM for the study area provided data sufficient for modeling groundwatersheds, groundwater flow paths, and the mapping of recharge and discharge. These features were the sole source of data used for interpolating the water-table surface and for estimating the depth to the base of the unconfined aquifer. While wetlands with specific hydrologic codes and hydric soil polygons were selected in an effort to identify the most persistent water-table elevations, it is possible that a portion of the wetlands included in the dataset may be perched or disconnected from the water table and do not have water elevations that truly reflect the permanent water-table surface.
Several areas where the gradient of the water table exceeds 0.05 were identified. Springs that were sampled for this study, along with those included in the Minnesota Spring Inventory [57] tend to occur at the lower margin of the slopes. These springs likely reflect the terrain where the water table intercepts the surface or focused flow related to the occurrence of less permeable deposits [19]. The springs could also be caused by preferential flow through high-permeability material [3].

Groundwater Flow Paths and Field Measurements
Measurements of pH varied widely in the study area and did not correlate with groundwater flow path lengths. A variety of factors influence surface water pH, including aquatic plant productivity, biological activity, and subsurface material [21,22,27]. As aquatic plant productivity and surface water pH in ponds is coupled with diurnal processes [22], it is possible that correlation may have been observed if pH measurements could have been taken at the same time of day at each site, preferably in the early morning hours. This would have significantly increased field sampling days and expense for this study, however.
Specific conductance correlated with groundwater flow path lengths at the 95% confidence limit and was a useful approach in this case. These results suggest lakes and wetlands lying on distal portions of groundwater flow paths are more mineralized and have water budgets that are more strongly influenced by groundwater; those lying on proximal portions of a groundwater flow paths are less mineralized and have water budgets that are not strongly influenced by groundwater.
A few exceptions in the model need to be highlighted. Lakes that lie at the beginning of flow path lines, including Augusta Lake and Horn Lake, had large specific conductance relative to their path line position. In addition, Nicolet Lake, which had the highest specific conductance measurement in the study area, lies toward the middle of a modeled groundwater flow path. These specific conductance measurements suggest that these lakes are at the end or middle of local flow paths, which is contradictory to the modeled flow paths. It is possible that Nicolet Lake is being fed by significant seeps or springs. The area of its groundwatershed is approximately nine times larger than the area of its surface watershed. The lake also lies at the base of the steeply sloping zone of the water table and a groundwater spring has been mapped just south of the lake. Unnamed ponds in the northwest portion of the study area had the smallest specific conductance measurements but lie toward the middle of modeled flow paths. The specific conductance measurements suggest these ponds are at the beginning of local flow paths but again contradict the modeled flow paths. It would be useful to sample ponds farther upgradient to determine if they have even lower specific conductance.
The distance groundwater travels in the subsurface is only one of many factors that can influence water quality, including specific conductance. For example, the residence time is also important to consider because the longer water remains in the subsurface the more likely its dissolved solids concentration will increase [58]. The composition of the subsurface materials and the sequence at which the groundwater encounters these materials are also important to consider as certain soil, rock, or sediment types may have minerals that are more prone to dissolution than others [58]. For example, more abundant shale fragments in the St. Louis sublobe likely influence mineralization differently than the otherwise calcareous drift. While the mineralization of water is a natural occurrence in the subsurface, it can be enhanced by human activities. For example, the use of salt to deice roads [59] or other human activities such as the extraction of natural resources [60] can increase the specific conductance of surface water bodies.

Watersheds and Field Measurements
Neither pH nor specific conductance were found to correlate with the ratio calculated from the surface watershed and groundwatershed areas (SW:GW). Although more than half the sampled lakes and wetland had surface watersheds and groundwatersheds that differed significantly in size, the surface watershed to groundwatershed ratio (SW:GW) could not be used to predict pH and specific conductance results. While basin size has been shown to positively correlate with residence time [28], which is an important factor in mineralization [3] and pH change [28], other properties of the watershed can influence the water budgets of lakes and wetlands and their water quality. For example, the surface watershed contribution is subject to soil properties, precipitation amounts, and losses from evapotranspiration. The groundwatershed contribution is subject to subsurface properties and recharge and discharge rates [3,61]. It is still useful to delineate both surface watersheds and groundwatersheds, however, especially for management purposes. For example, there can be large differences between the two (e.g., Figure 12) and knowing the catchment area for both is useful when considering land use changes, groundwater development, and contaminant transport. It is likely that the areal extent of groundwatersheds differs between dry and wet conditions due to fluctuations of the water-table surface that result from changes in recharge and discharge [3]. This is also an important consideration for both modeling and management, especially when taking into account the implications of climate change.

Recharge-Discharge Map
Strong discharge areas in the study area appear to correspond to mapped groundwater springs and areas of the water table with steep gradients. It would be useful to explore other strong discharge areas in an attempt to identify additional springs. In addition, one wetland sampled east of Whipple Lake, believed to be an ombrotrophic bog due to its low pH and vegetation, mapped as a strong recharge area surrounded by discharge. This type of recharge-discharge pattern could also be further explored in an attempt to identify additional bogs within the Itasca Moraine and elsewhere.
For the recharge-discharge map, it was assumed that the Ksat value obtained from the SSURGO soils [38] was a reasonable estimate for hydraulic conductivity of the shallow glacial sediments and that the unconfined aquifer is laterally isotropic. As there are no observational or test data, the Ksat value obtained from the SSURGO soils, mapped at a scale of 1:20,000 and 1:24,000 [38], is the best source for the hydraulic conductivity of surficial sediments available in the study area. In addition, no Ksat values were available for lake sediments within the study area so they were assigned a Ksat value similar to that of wetlands in the study area. According to Wright [30], the organic sediments of lakes in the area can range from 8 to 12 m thick. The hydraulic conductivity assigned to the lakes is in range of what Mięsiak-Wójcik et al. [62] found in lake bottoms composed of fine organic sediment with low percentages of fine sand. The aquifer base elevation, which we calculated using an analytical model developed by Bresciani et al. [52], used the Ksat values [38] and recharge estimates [54]. Therefore, the aquifer base elevation was not based on an observed confining layer, for which data are unavailable, but instead on model results.
The computed recharge and discharge rates shown in the recharge-discharge map ( Figure 13) are generally larger than regional estimates [54]. Large gradients, which may generate the large modeled groundwater flux, especially on the west side of the watershed (Figure 13), exceed rainfall rates by as much two orders of magnitude. This may be partly explained by the redistribution of surface water by surface flow. More likely, however, large fluxes may result from excessively large hydraulic conductivity estimates based on soils ( Figure S4) that do not account for layer anisotropy, decreasing hydraulic conductivity with depth in glacial sediments, and/or aquifer saturated thickness values derived from a constant aquifer base elevation. This condition may occur due to perched water tables and greatly diminished hydraulic conductivity as a result of organic and fine sediment accumulation in wetlands and lakes. For example, Rosenberry [63] observed unsaturated conditions and large hydraulic flux (measured as high as 2.63 m/day) beneath Lake Belle Taine, about 30 km southeast. Due to its uncertainties, the recharge-discharge map ( Figure 13) provides a semiquantitative estimate of recharge and discharge rates in the study area, although the pattern of recharge and discharge across the groundwatershed is reliable.
Future work with this method could use a constant-thickness, hydrologically active zone for the aquifer rather than relying on a fixed aquifer base elevation across the entire study area. This approach would reduce the large residuals that occur in more elevated areas where the aquifer thickness is greatest. It is possible that the values estimated from the soils are too high and do not adequately represent the more compacted underlying sediments [64,65]. Deep borings and data supplied by slug or pump tests would be useful to refine the hydraulic conductivity and the aquifer thickness estimates, but would obviously increase the research cost significantly.
Sensitivity tests indicate that the 40-m cell size and the raw (unfilled) water-table surface provided the best match to the discharge measurements from the two outlets (Table 2) and the overall water budget for the watersheds. Although this is coarser than the DEM used to determine the surface watersheds, rasters with a coarser resolution are likely more appropriate for models of the water-table surface. This is based on the assumption that the water-table surface is smoother and more subdued than the surface topography [11,28]. While filling the water-table surface was appropriate for delineating groundwatersheds, it increased net residual discharge for the study area in the groundwater volume balance residual raster. This is because filling the surface sinks results in small increases in the heads of various cells [47]. The difference between the discharge of the filled and raw (unfilled) water-table surface, integrated over the entire groundwatershed, was as little as 0.02 m 3 /s for the 40-m cell-size trial and as much as 1.3 m 3 /s for the 30-m cell-size trial. As DEM surfaces are filled to ensure a connected drainage structure and therefore enable delineation of watersheds and groundwater flow path lines, using a filled raster for the production of the recharge-discharge map is probably unnecessary. In addition, these areas could represent actual groundwater sinks and not model artifacts, thus warranting further investigation.

Specific Conductance, Flow Paths, and Mapped Glacial Sediments
When comparing the flow paths and average specific conductance measurements to the mapped glacial deposits, measurements of <100 µS/cm were recorded in lakes and wetlands located on the more elevated areas of the moraine, which are composed of glacial till. The bog near Whipple Lake is an exception and lies within an area of generally greater conductance. Specific conductance measurements >100 µS/cm were measured in the lakes and wetlands situated on the glacial outwash or remnant tunnel valley features within the study area. Similarly, in the Young Glacial Area of northern Poland, Jaworksa-Szulc [13] found that recharge or losing lakes with low total dissolved solids (TDS) were commonly identified on till while discharge or gaining lakes with high TDS were commonly identified within glacial tunnel valleys. These results show the importance of hydraulic conductivity and position within the groundwater flow system for characterizing surface water and groundwater interaction [3].

Specific Conductance Variability Compared to the Missouri Coteau and Young Glacial Area
Lakes and wetlands of the Cottonwood Lake area of the Missouri Coteau, North Dakota, and the Young Glacial Area of northern Poland have been studied in an effort to determine their interaction with groundwater. Although these areas and the Itasca Moraine consist of similar glacial terrain, there are notable differences in climate and characteristics of surface water mineralization.
The "prairie potholes" in the Cottonwood Lake area have been studied extensively. In particular, many specific conductance measurements were collected for 33 years from 1979 to 2012 [66]. Six Cottonwood potholes (P1, P8, T3, T4, T5, T8) exhibiting a range of geochemical characteristics were selected for comparison to the Itasca site ( Figure 14). Ponds P1 and P8 are described as discharge ponds, T3 and T4 as flow-through ponds, and T5 and T8 as recharge ponds [67]. In contrast to the Itasca Moraine, evapotranspiration exceeds precipitation in the Cottonwood Lake area [68], and wetland waters, derived from precipitation and snowmelt, on average, are more mineralized than those found along the Itasca Moraine. The higher evapotranspiration rate in this area influences specific conductance as it causes salts to accumulate in the ponds. The ponds in this area are believed to lose little water to groundwater recharge, and the groundwater input to the ponds is believed to be small [69,70]. In cases where the ponds do interact with groundwater, dissolved minerals can be flushed from the ponds through groundwater recharge and large amounts of minerals can be supplied to the ponds through groundwater seepage [71]. The two discharge ponds selected, P1 and P8, had an average specific conductance of 1900 µS/cm. The two flow-through ponds, T3 and T4, had an average of 1400 µS/cm while the two recharge ponds, T8 and T5, had an average of 190 µS/cm.
The glacial lakes of the Young Glacial Area of northern Poland have total dissolved solids (TDS) measurement data from 2010 to 2012 by Jaworska-Szulc [13]. Of these lakes, 15 were described as being discharge lakes, three as flow-through lakes, and seven as recharge lakes. For comparison, TDS measurements by [13] were roughly converted to specific conductance by dividing TDS values by 0.6 [72] (Figure 14). Precipitation in the Young Glacial Area nearly equals the evapotranspiration rate [13] and the lakes are, on average, less mineralized than those found in the Cottonwood Lake area. The 15 discharge lakes (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15) had an average calculated specific conductance of 460 µS/cm. The three flow-through lakes (26,34,35) had an average of 180 µS/cm and the seven recharge lakes (16)(17)(18)(19)(21)(22)(23) had an average of 110 µS/cm. Figure 14. Specific conductance measurements for the Cottonwood Lake area (North Dakota) [12], Young Glacial Area (northern Poland) [13], and Itasca Moraine (this study) using averaged values for the types of ponds and lakes sampled. In the plots, the x symbol represents the mean, the middle bar the median, the upper box the third quartile, the lower box the second quartile, the upper whisker the fourth quartile and the maximum, and the lower whisker the first quartile and the minimum.
The gaining lakes at the Poland site are the most mineralized, the flow-through lakes have intermediate mineralization, and the losing lakes have the least mineralization [13]. This corresponds to the general trend at the Cottonwood Lake and Itasca sites, although the magnitude of mineralization is strikingly different among the three sites ( Figure 14).
In this comparison, it is important to emphasize that different factors influence mineralization. The Cottonwood Lake area has evapotranspiration rates that exceed precipitation rates by approximately 30 cm/year [68], whereas these rates are approximately equal for the Young Glacial Area [13], and precipitation exceeds evapotranspiration by up to 5 cm/year in the Itasca region [73]. It is possible the that lower evapotranspiration rates and the increased interaction with groundwater occurring at the Itasca Moraine and the Young Glacial Area result in a flushing effect that has kept dissolved minerals from becoming concentrated in the lakes and wetlands. The Cottonwood Lake area ponds are described as having less interaction with groundwater relative to their water budget, which is one mode for mineral import and export. This site also experiences greater evapotranspiration which can mineralize pond water through concentration processes [69,70]. When water levels are high, the potholes fill with water and can spill into adjacent basins, which also exports minerals [74]. In addition to variations in the composition of sediment, Heagle et al. [75] suggested that past dry periods in Saskatchewan's prairie pothole region caused a large mass of sulfate to accumulate in the wetland's sediments and that the sediments release ions to the surface water during wet periods, thereby maintaining salinity. Figure 14. Specific conductance measurements for the Cottonwood Lake area (North Dakota) [12], Young Glacial Area (northern Poland) [13], and Itasca Moraine (this study) using averaged values for the types of ponds and lakes sampled. In the plots, the x symbol represents the mean, the middle bar the median, the upper box the third quartile, the lower box the second quartile, the upper whisker the fourth quartile and the maximum, and the lower whisker the first quartile and the minimum.
The gaining lakes at the Poland site are the most mineralized, the flow-through lakes have intermediate mineralization, and the losing lakes have the least mineralization [13]. This corresponds to the general trend at the Cottonwood Lake and Itasca sites, although the magnitude of mineralization is strikingly different among the three sites ( Figure 14).
In this comparison, it is important to emphasize that different factors influence mineralization. The Cottonwood Lake area has evapotranspiration rates that exceed precipitation rates by approximately 30 cm/year [68], whereas these rates are approximately equal for the Young Glacial Area [13], and precipitation exceeds evapotranspiration by up to 5 cm/year in the Itasca region [73]. It is possible the that lower evapotranspiration rates and the increased interaction with groundwater occurring at the Itasca Moraine and the Young Glacial Area result in a flushing effect that has kept dissolved minerals from becoming concentrated in the lakes and wetlands. The Cottonwood Lake area ponds are described as having less interaction with groundwater relative to their water budget, which is one mode for mineral import and export. This site also experiences greater evapotranspiration which can mineralize pond water through concentration processes [69,70]. When water levels are high, the potholes fill with water and can spill into adjacent basins, which also exports minerals [74]. In addition to variations in the composition of sediment, Heagle et al. [75] suggested that past dry periods in Saskatchewan's prairie pothole region caused a large mass of sulfate to accumulate in the wetland's sediments and that the sediments release ions to the surface water during wet periods, thereby maintaining salinity. These are factors that may contribute to the higher mineralization found in the prairie potholes.
It is also important to consider land use when comparing sites. The Itasca location is largely undisturbed forest while the Cottonwood Lake area is primarily native prairie [66]. The Young Glacial Area is a mix of forested and agricultural land [13]. Depending on the type of land use or cover, the runoff from the land surface can increase the specific conductance in lakes and wetlands. Climate change may also influence differences in mineralization among the sites. Elk Lake in Itasca State Park has been the focus of many paleoclimate studies as it lies within the transient North American forest-prairie ecotone. Elk Lake is unusually deep (>30 m) and has well-developed annual varves that have provided key information on past climate variability [76][77][78][79]. The results of studies on these varves suggest that the Itasca area is sensitive to climate change and was a prairie ecosystem during the warm and dry middle Holocene [76]. Although the lakes and wetlands at Itasca are least mineralized among the three sites, it is possible that during the Holocene these lakes and wetlands were similar hydrogeochemically to current conditions in the Cottonwood Lake area.

Conclusions
Application of this GIS-based groundwater modeling technique to the Itasca Moraine provided valuable water budget information and would likely produce good results in other areas with abundant streams, lakes, wetlands, and hydric soils. Having subsurface geological data would strengthen model outcomes. Flow path length and specific conductance showed a significant correlation, suggesting that this relatively simple and inexpensive technique that couples GIS modeling and field measurements can help characterize the contribution of groundwater in wetlands.
In contrast to pathline length, the ratio of watershed areas was not shown to be significantly correlated to either pH or specific conductance, knowing the catchment for both surface water and groundwater is still imperative to make informed management decisions. This is because the two catchments do not always correspond [5] and changes to either of the catchments, in the form of land use or groundwater development, for example, can directly affect the quantity and quality of water within a lake or wetland. Mapping recharge and discharge also can be of use to management for identifying groundwater springs and sensitive wetlands. In addition, knowing the distribution of strong recharge and discharge sites would be useful when making decisions that affect land use within a watershed.
The Cottonwood Lake prairie site in North Dakota (U.S.) and the Young Glacial Area (northern Poland) also show a trend of increased specific conductance in lakes and wetlands with increased groundwater interaction. Results for the Itasca Moraine suggest that this GIS modeling technique would function similarly in those areas, even with contrasting precipitation-evapotranspiration ratios and different magnitudes of surface water mineralization. In addition, the prevalent land use or land cover of the watershed and the nature of sediments in the wetlands must be considered in the GIS analysis as they are shown to influence specific conductance measurements of surface water.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/w13070983/s1, Figure S1: Three-meter resolution digital elevation model of the study area, Figure S2: Profile showing the effect of adding additional water-table control features in an area of complex topography, Figure S3: Wetland and lake polygons used to interpolate the water-table surface, Figure S4: Saturated hydraulic conductivity Ksat for the study area based on SSURGO soil polygons and weighted average soil series profiles, Figure S5: Depth to the interpolated water-table surface for the study area, Figure S6: Scatter plot of flow path lengths and pH measurements for sampled lakes and wetlands, Figure S7: Scatter plot of the watershed area ratio and pH measurements for sampled lakes and wetlands, Figure S8: Scatter plot of the watershed area ratio and specific conductance measurements for sampled lakes and wetlands.    [52] developed a 2D topographic profile model that estimates the recharge-hydraulic conductivity ratio of a region based on measurements of elevation and distance to the nearest constant head boundary ( Figure A1). We applied the model to estimate the R/K ratio for the lakes and wetlands in the Itasca Moraine in Northern Minnesota. Based on USGS estimates for groundwater recharge [54] and SSURGO data [38] to calculate average hydraulic conductivity the expected R/K ratio for the study area ranges from 3 × 10 −5 to 2 × 10 −4 m/day. After adjusting the maximum depth of surfacegroundwater interaction, the model output R/K ratios were achieved within the expected range. Additional stratigraphic data are needed to improve the estimated depth limit to interactions between surface and groundwater.
The Bresciani et al. model [52], which relies on the Dupuit assumption, relates recharge and hydraulic conductivity to where the water table meets the surface along a horizontal profile. The five variables shown in Figure A1 were selected and tested by [52] using MODFLOW-NWT [50,80,81] to yield an R/K value, which resulted in the water table reaching the surface only at measured surface water features [52]. The resulting numerical solution is a modification of Hooghoudt's drainage equation [53]. The model includes three simplifying assumptions: groundwater flow is unconfined and two-dimensional, the aquifer is homogeneous, and groundwater flows within the cross-section between the two constant head boundaries. Used in this context, the constant head boundaries of the profile are open bodies of water that are assumed to have a fixed hydraulic head.

B. Theory and Methods
Several equations and coupled relationships comprise the Bresciani et al. model [52]. Equation (A1) relates the geometry of the topographic profile and water-table elevation ( Figure A1) to the average ratio of recharge rate to hydraulic conductivity (variables given in Table A1) [52]: where the equivalent depth De is less than D, which accounts for the vertical flow resistance below the level of the drain. De is given by [52]: where ro (L) is the radius of the drains and 2πD

B. Theory and Methods
Several equations and coupled relationships comprise the Bresciani et al. model [52]. Equation (A1) relates the geometry of the topographic profile and water-table elevation ( Figure A1) to the average ratio of recharge rate to hydraulic conductivity (variables given in Table A1) [52]: where the equivalent depth D e is less than D, which accounts for the vertical flow resistance below the level of the drain. D e is given by [52]: where r o (L) is the radius of the drains and u = 2πD L (A3) Equation (A2) adjusts D to allow for resistance to vertical groundwater flow below the elevation of the reference drain.
This converges rapidly when u > 0.5. If u < 0.5, van der Molen and Wesseling [82] recommend using an approximation from Dagan's formula [83]:  (A4)) is a function of u, which is a function of the shape of the drain (Equation (A5)). These equations and variables, described in greater detail by [52], were programmed into an Excel spreadsheet for ease of repeated calculations.
Several steps were used to apply the method. First, constant head boundaries were selected within the Itasca State Park area. Historical satellite imagery from 1991 to 2013 from Google Earth was examined to find surface water features with visually consistent water levels through wet and dry years. Elk Lake, Lake Itasca, and Basswood Creek were selected as constant head boundaries because they exhibit minimal fluctuations (less than 0.5 m) in total hydraulic head seasonally and from year to year (Table A2). Next, the topographic profiles for the transects noted in Table A2 were derived from the LiDAR digital elevation model using QGIS and SAGA terrain analysis [84]. Areas where surface water features intersected with each profile were identified visually with the satellite imagery and Minnesota wetland delineation spatial data [35,85]. Wetlands within 500 m, or approximately 5% of the total profile length, on either side of each profile were projected onto the profile to add more data points where few wetlands directly intersected the profile. Each profile was used to measure L, x, z t , and H. Variable D was initially approximated using a map of the estimated depth to bedrock [86], under the assumption that no significant confining layers are present between the surface and bedrock. Variables x and z t were input for each surface water feature along the profile into the ratio calculator spreadsheet to yield separate R/K values for each surface water feature.

C. Results and Discussion
Using the SSURGO Web Soil Survey interactive tool [38] to select an area of interest along each of the three measured profiles, the estimated average hydraulic conductivity along the three measured profiles ranged from 2.3 to 9.6 m/day. Due to lack of data for depths greater than 1.5 m (5 feet), this calculation assumes that the hydraulic conductivity near the surface is similar to subsurface hydraulic conductivity for the entire depth of D. These hydraulic conductivity and recharge estimates yield an expected R/K ratio range from 4 × 10 −5 to 2 × 10 −4 (Table A3). Table A3. Expected values for average hydraulic conductivity (K), recharge (R), and R/K ratios based on weighted average K [38] and recharge rates from [54].

Profile
Average K (m/day) R (m/day) Expected R/K Initial calculations yielded average R/K values of 2.6 × 10 −4 , 3.2 × 10 −4 , and 3.6 × 10 −4 for profiles 1, 2, and 3, respectively (Table A4), which are all larger than the expected R/K range. The calculated R/K value for each profile is roughly an order of magnitude larger than the expected R/K value. By adjusting the value for depth to the shallowest confining layer beneath the profile (D), the calculated R/K ratios were achieved within the expected range.
This application of the model revealed its sensitivity to variable D, the depth to a horizontal confining layer (e.g., Figure A2). Initially, variable D was assumed to be equal to depth to bedrock for all profiles, or approximately 200 m. However, the expected R/K values may suggest much shallower depths to a confining layer at 5, 70, and 30 m for Profiles 1, 2, and 3, respectively (Table A4). Glacial till can be heterogeneous, and soil series such as Haslie, Seelyeville, and Cathro, which are present in the region, may have a sufficiently low hydraulic conductivity to act as a confining layer.  This application of the model revealed its sensitivity to variable D, the depth to a horizontal confining layer (e.g., Figure A2). Initially, variable D was assumed to be equal to depth to bedrock for all profiles, or approximately 200 m. However, the expected R/K values may suggest much shallower depths to a confining layer at 5, 70, and 30 m for Profiles 1, 2, and 3, respectively (Table A4). Glacial till can be heterogeneous, and soil series such as Haslie, Seelyeville, and Cathro, which are present in the region, may have a sufficiently low hydraulic conductivity to act as a confining layer. Figure A2. The nonlinear relationship between depth to the horizontal confining layer (D) and the recharge/hydraulic conductivity ratio (R/K). Note that the effect on R/K is diminished as D becomes greater than about 200 m. D values varied across profiles, but so did expected hydraulic conductivity values (Table A4). This variation implies either the depth to which groundwater interacts (D) is not constant within the profile, or the measured hydraulic conductivity (K) was too detailed for the dimensions of this model. Applications of the model should consider incorporating a more general K value, and account for uncertainty caused by variations in D. Figure A2. The nonlinear relationship between depth to the horizontal confining layer (D) and the recharge/hydraulic conductivity ratio (R/K). Note that the effect on R/K is diminished as D becomes greater than about 200 m. D values varied across profiles, but so did expected hydraulic conductivity values (Table A4). This variation implies either the depth to which groundwater interacts (D) is not constant within the profile, or the measured hydraulic conductivity (K) was too detailed for the dimensions of this model. Applications of the model should consider incorporating a more general K value, and account for uncertainty caused by variations in D.
This application revealed the Bresciani et al. model's sensitivity to the depth to a confining layer (D) [52]. Initially, D was assumed to be equal to depth to bedrock, due to limited geologic log data from boreholes in the Itasca State Park region. However, the results indicate that the depth limit of surface-groundwater interaction is much shallower, possibly roughly correlating to the depth of Elk Lake or Lake Itasca, which was applied to the groundwater recharge-discharge model described in the main text.