Interactive comment on “ Hydrogeological controls on spatial patterns of groundwater discharge in peatlands ”

“Hydrogeological controls on spatial patterns of groundwater discharge in peatlands” by Hare et al. describes a multifaceted approach to understanding two primary mechanisms for groundwater discharge in peatlands: matrix seepage and preferential flow path seepage. The research is primarily based in temperature observations (DTS, vertical temperature profilers, and infrared imaging), but also incorporates regional well data, ground penetrating radar, and coring information. The authors conclude that peatlands operate as part of the regional groundwater system, and that underlying basin curvature and peat thickness likely control patterns of discharge.


Introduction
Peatlands develop in response to physical, biological, and chemical processes and feedbacks.Groundwater discharge to surface water is one of the most important physical controls on peatlands stability (Siegel et al., 1995;Watters and Stanley, 2007); yet the underlying physical hydrogeologic framework governing the development of surface seepage distribution in these systems is not well understood.Preferential flow paths, hydraulic conductivity (K) anisotropy, and geologic heterogeneities likely influence the surface expression of discharge zones (Chason and Siegel, 1986;Drexler et al., 1999;Smart et al., 2012).However, these variables have been difficult to constrain due to the spatial resolution of traditional localized groundwater wetland methods (wells, boreholes, surface point measurements, etc.) and their impact on fragile flow paths.The underlying hydrologic engine of these Published by Copernicus Publications on behalf of the European Geosciences Union.D. K. Hare et al.: Hydrogeological controls on spatial patterns of groundwater discharge wetlands have shown to be difficult to discern in large-scale systems.
Thermal dynamics of ground and surface waters also govern critical wetland functions and can be assessed in multiple ways.Surface water thermal stability, for example, is a popular research focus in ecohydrology, as this process is important for aquatic species that rely on the low variance of groundwater temperature to buffer themselves from heat extremes and regulate their metabolism (Caissie, 2006;Deitchman and Loheide II, 2012).Temperature also controls chemical processes in ecosystem respiration, which in turn controls carbon processing and nutrient retention (Boulton et al., 1998;Davidson and Janssens, 2006;Demars et al., 2011;Lafleur et al., 2005), biodiversity (Parish et al., 2008), as well as overall species health (Verberk et al., 2011).Upwelling zones are linked to increased biogeochemical cycling (Sebestyen and Schneider, 2001) and also maintain species richness through the "edge effect" -overlap between the thermally and chemically stable groundwater ecotone and the higher oxygen environment within the main stream channel (Brunke and Gonser, 1997;Cirkel et al., 2010).An increase in wetland temperature has been shown to stimulate methane production (McKenzie et al., 2007) as well.The underlying drivers of the thermal regime of a wetland system can be caused through varying driving processes, and are important to the ecosystem services provided in the peatland.
Widespread drainage of peatlands has caused wetland degradation and loss of ecosystem services.Anthropogenic modifications such as ditching and filling create discontinuity between surface water and groundwater systems, with impacts on wetland function (van Loon et al., 2009).In some parts of the world, wetland restoration is attempting to address these historical impacts.Within the United Kingdom, for example, efforts to return natural water table levels by filling drainage ditches in peat mining areas have led to disagreements as to the cost-benefit of these specific restoration designs (Grand-Clement et al., 2013).In New England (United States), where thousands of acres of historical peatlands were converted to commercial cranberry farming in the late 1800s (Garrison and Fitzgerald, 2005), wetland restoration is similarly attempting to regain natural water table levels (Price et al., 2003).An incomplete understanding of the underlying hydrology and thermal regime can limit the effectiveness of such efforts.
In this research, we explore the spatial distribution of groundwater seepage through a kettle-hole peatland from the analysis of basin structure and hydraulic properties of the peatland matrix.To assist in wetland restoration design at the study site, we focus on understanding the natural processes that promote the hydrologic inputs for aquatic habitat formation and maintenance.The goals of this study are to (1) identify groundwater discharge locations and their hydrogeologic controls, (2) determine temperature dynamics of the groundwater discharge locations, and (3) evaluate the development of these seepage patterns.Through this work, in-sight is gained into how the hydrologic driving mechanisms of peat-based wetlands can support the restoration of sustainable ecosystems (e.g., process-based design) (Beechie et al., 2010;Dahl et al., 2007).

Site description
The site "Tidmarsh Farms" is comprised of three cranberry farms, and the two largest farms are separated by Beaver Dam Road (Fig. 1).The area surrounding this peatland site is characterized by outwash, kame deltas, and ground moraines that show evidence of collapse features and deformation (Larson, 1982;Stone et al., 2011).These ice collapse features are typical of environments proximal to ice contact zones and can result in the formation of kettle holes, of which there is extensive evidence throughout the surrounding region.All three of the site's cranberry farms were built on kettle-hole peatlands between the late 1800s and early 1900s.The cranberry peatlands on Tidmarsh East were taken out of production in 2010 and another was taken out of production in 2015.This study concerns work on Tidmarsh East (the site), which is 2.5 km 2 (Fig. 1b).
A layer of sand 0.3-1.5 m thick overlays native surficial soils, which was laid down as a part of a normal cranberry farming practices until the site retirement in 2010.This practice maintained a very low gradient across the site with a slight decline to the north, with minimal microtopography.
During this research, conducted in collaboration between the Living Observatory and University of Massachusetts, a restoration project involving the private landowners, governmental agencies, and nongovernmental organizations was in assessment and design phases.Project planners were specifically interested in the location of groundwater discharge across the site to help design the placement of reconstructed stream channels.In addition, the restoration design team sought to better understand the location of subsurface peat deposits, underlying site hydrology, and potential future thermal regimes when considering potential restoration activities.As of 2017, the site has undergone both passive and active restoration to encourage an accelerated ecological recovery based on the conclusions of this work.In the following sections, we document our methods and findings specific to the spatial distribution of the groundwater discharge at Tidmarsh East and the implications for restoration design.

Site hydrology
The farm is a part of the small 5 km 2 surficial Beaver Dam Brook Watershed, but is also a discharge location of the 360 km 2 Plymouth-Carver-Kingston-Duxbury (PCKD) aquifer, and thus the groundwater flow paths contribute from a much larger hydrologic system than the surficial watershed (Fig. 1a).The PCKD aquifer below the site is characterized by glacial outwash sands (Masterson et al., 2009).Surface water enters the site from four surface water bodies south Hydrol.Earth Syst.Sci., 21, 6031-6048, 2017 www.hydrol-earth-syst-sci.net/21/6031/2017/ Figure 1.(a) Site map of the Tidmarsh Farms regional peatland showing the study area and watershed boundary; Plymouth County, Massachusetts, and PCKD USGS wells used for regional groundwater isotopic data (Table 1).(b) Detail of the Tidmarsh Farms study site showing the major waterways and flow direction in blue, site groundwater wells, isotopic sample locations, and GPR transects.Beaver Dam Brook flows north into Plymouth Bay.
of the site (Fresh Pond, Little Island Pond, the Arm Wetland, and Beaver Dam Pond headwaters) and drains northward into Beaver Dam Brook, an approximately 2 km reach, before discharging in Bartlett Pond and then directly into Plymouth Bay (Fig. 1b).
To facilitate drainage and irrigation, lateral and perimeter drainage ditches exist throughout farmed areas.Parallel drainage ditches are located approximately every 18-35 m throughout the entire site, and are approximately 1 m wide and 0.5 m deep.The western agricultural cells have drainage ditches oriented east-west (Cell 3 and 4), and in the eastern cells (Cell 6 and 7) most drainage ditches are oriented north-south.When the study was conducted, the site was still predominately covered in low-lying cranberry vegetation, as well as a variety of sedges and cattails mostly adjacent to the central stream bank and marginal drainage ditches.
Flashboards in the dam creating the Beaver Dam Pond impoundment were permanently removed by the landowners in the fall of 2010, after which the southern side of the farm was allowed to return to a natural wetland state (Fig. 1b).Data collection conducted for this research spanned 2012-2014, beginning 2 years after farming ceased, and prior to any active wetland restoration activity.
The site is located within the discharge zone of the large PCKD aquifer, and thus short-term, drastic temporal shifts are not expected in the hydrogeology or the processes described herein.We expect that our observations from the study conducted over this 2-year period to be representative of present-day conditions.The primary source of recharge to the PCKD aquifer is through precipitation which rapidly infiltrates outwash plain deposits (Wareham and Carver Pitted Plains) (Masterson, 2009), and thus changes in the water table elevation can be expected.Topographic changes to the base level due to isostatic rebound and sea level rise may also contribute to water table elevation changes (Oakley and Boothroyd, 2012).The regional aquifer may be sensitive to long-term climatic changes (Shuman et al., 2001;Newby et al., 2009); however, this question is outside the scope of this study.

Methods
Seepage patterns within peatlands have been difficult to constrain due to large site areas and complex, dynamic substrates.At Tidmarsh Farms, we use multiple remote-sensing and direct-contact methods in this environment to connect data from different scales into a process-based understanding of peatland groundwater seepage.Ground-penetrating radar (GPR) is used to evaluate the subsurface structure of the peatland basin(s), and multiple thermal methods are used to locate and analyze surficial groundwater seepage patterns.Stable water isotopes are used to describe the dominant water sources supplying the seepage.
Traditional hydrogeologic methods were also implemented including well transects, seepage meters, and differential discharge gauging along Beaver Dam Brook.

Resolving subsurface structure
GPR has been successfully used to characterize peatlands' physical structure and stratigraphy due to the strong contrast between peat and the underlying aquifer geophysical properties (e.g., water content) (e.g., Comas et al., 2005;Holden, 2004;Kettridge et al., 2008;Lowry et al., 2009;Slater and Reeve, 2002).The GPR method relies on the transmission of electromagnetic (EM) waves through the subsurface then records the time and amplitude of the returning signal (reflection) to image changes in the EM properties between subsurface materials (Knight, 2001;Lowry et al., 2009).In August 2012 we collected common-offset reflection data using Malå ProEx 100 and 50 MHz antennas with a transmitter-receiver separation of 1 and 2 m, respectively.Here, we only use the 100 MHz data to generate interpolated maps of peat thickness, as those data provide better resolution of the peat-sand interface given the EM properties of the peat matrix and the depth of the structure of interest, which was 0-15 m for this study site.A total of 19 GPR line surveys were completed; all surveys used 0.3 m trace spacing and ranged from 100 to 1000 m in total length (Fig. 1b).The vertical resolution of the survey was 0.9 m, based on the theory that layers can only be resolved if their thickness is greater than one-quarter of a wavelength.
We applied a 150 MHz high-cut filter to remove the highfrequency noise, and then a 100 ns automatic gain control to compensate for signal loss with depth, and distinguish deeper reflections by averaging over the time window applied and adjusting the central signal strength with respect to that result.No topographic adjustments were made, as there is negligible topographic variation both along the surveys and between the surveys.The peat thickness was determined in each of the radargrams.Three characteristic radargrams are shown in Fig. 2.
We constrain the EM signal velocity through the peat for the GPR data analysis, and describe the peat's structure with depth by collecting nine sediment cores (Fig. 1b) with a vibracorer (3) and hand corer (6).Analysis of the cores demonstrated that the layered reflections observed in the radargrams were due to variation in the degree of humification.We determined an average EM velocity of 0.036 m ns −1 (range = 0.030-0.040m ns −1 ) through the peat for the five full length cores that extended to the peat-sand interface.This velocity range is consistent with other peatland GPR studies (0.033-0.040 m ns −1 ) (Parsekian et al., 2012).Using these data, a 3-D interpolation of the peat-sand interface was created using kriging to estimate the subsurface peat basin structure (Fig. 2).The second derivative of the maximum slope (profile curvature) was calculated from the interpolated surface to identify changes in basal slope of the peat-sand interface, and is shown in Fig. 2.

Identifying locations of groundwater discharge to surface water using temperature
Heat can be used as a tracer to identify upwelling groundwater, as air temperature oscillations on diurnal and annual timescales strongly influence surface waters, while deep (e.g., greater than approximately 10 m) groundwater temperatures remain relatively constant through time (Anderson et al., 2005;Constantz, 2008).Local, shallow flow paths can be more sensitive to climatic and seasonal changes in evaporation and precipitation (Fraser et al., 2001;Kurylyk et al., 2014b;Menberg et al., 2014;Reeve et al., 2006) and may not contribute to the thermal stability of aquatic systems to the same extent as deep (> 10 m) regional aquifers.This noted, during the thermal study periods, groundwater temperatures range from 10-11 • C in on-site wells below the peat.

Fiber-optic distributed temperature sensing
Raman spectra fiber-optic distributed temperature sensing (FO-DTS) is used for spatially extensive heat tracing in aquatic systems.Tyler et al. (2009) provides a thorough review of the details of the technology and calibration.DTS temperature data were collected with a SensorTran Gemini HT control unit in dual-ended mode using an AFL telecommunications umbilical fiber-optic cable.This FO-DTS unit allows for 1 m spatial accuracy at 0.1 • C precision over ∼ 15 min integration times.Each FO-DTS deployment was operated for a minimum of 5 days to ensure multiple sufficiently strong diurnal oscillations were captured.Calibration coils that are 50 m long were maintained at a constant temperature with an ice-water slush bath and/or ambient bath and were compared to an independent Onset HOBO Water Temperature Pro v2 Data Logger (U22-001) (±0.2 • C accuracy).
In July and August of 2013 four FO-DTS deployments were installed, one within the drainage ditches of eastern peatland cells, and three within the western cells.We capitalize on the modified structure of the agricultural peatland surface, particularly the relatively evenly spaced drainage ditches, to thermally sample surface water in a distributed way which is not possible in more natural systems (e.g., Lowry et al., 2007).The deployment sites were chosen based on previous infrared surveys (27 November 2012, discussed in Sect.2.2.2), interviews with the farmer, and feasibility of installation.Each deployment ranged from 1000 to 2500 m in length.Macrophyte growth was cleared during installation and continuously monitored through each deployment.
The arithmetic mean and standard deviation were calculated for each ∼ 5-day time series of FO-DTS data at every meter along the fiber-optic cable to identify locations of groundwater seepage.These results can indicate the location and relative magnitude and permanence of groundwater discharge, which is not possible with other methods, such as thermal infrared (TIR) or temperature probes (Briggs et al., 2012;Hare et al., 2015;Sebok et al., 2013;Selker et al., 2006).

Infrared surveys
TIR cameras sense and quantify surface infrared (heat) radiation and are increasingly being used to evaluate aquatic systems efficiently on large scales (Chen et al., 2009;Deitchman and Loheide, 2009;Dugdale et al., 2016;Handcock et al., 2012;Hare et al., 2015), particularly at large sites, or sites where in situ measurements are not possible.The hand-held TIR survey was conducted to both expand the thermal survey and to compare this method to the FO-DTS data.We used a high-resolution forward-looking infrared camera (T640BX model FLIR, FLIR Systems, Inc.) with GPS and compass capabilities.The TIR method allowed for efficient spatial coverage and allowed us to obtain thermal data unreachable with FO-DTS (Hare et al., 2015).
At Tidmarsh Farms East three TIR surveys were completed: one spanning 30-31 July 2013, one on 21 March 2014, and one reconnaissance survey on 27 November 2012.The July survey was used to make comparisons to the FO-DTS data as it was taken during the same time period; the March survey was used to compare seasonal variability in seepage patterns.Surveys were conducted in the morning and evening to minimize reflection interference, and all temperature collection practices and considerations for this site are described in detail in Hare et al. (2015).To create a spatial site map comprised of all TIR images, a single temperature (color-contoured pixel) from an aquatic point of interest was selected, and used to color an icon on the map.This allowed for georeferenced TIR data to be used quantitatively to evaluate seepage patterns by location.The relative magnitude of the seepage rate is estimated based on how similar the observed temperature is to the regional groundwater temperatures.

1-D vertical temperature profiles
The depth to which the surface diurnal temperature signal penetrates saturated near-surface sediments depends on the period of the signal, the fluid flow velocity and direction, and the physical properties of the fluid-saturated sediment (Goto et al., 2005;Hatch et al., 2006;Irvine et al., 2017;Stallman, 1965).With depth, the diurnal heat signal variation decreases in amplitude and its shifts forward in time.Much of the heat transport not explained by pure conduction is attributable to advective fluxes, which can be solved for from thermal time series at multiple depths using simple analytical solutions to the 1-D heat transport equation with specified boundary conditions (Hatch et al., 2006;Rau et al., 2014;Schmidt et al., 2007;Silliman et al., 1995;Stallman, 1965).
We analyzed four 1-D vertical temperature profiles to understand the vertical subsurface fluid flux patterns at the site.Maxim iButton temperature loggers (0.0625 • C resolution; 1 • C accuracy) were attached to cavities drilled into a wooden dowel and placed into the ground such that the logger locations were −2.5, −5.0, −10.0, and −25.0 cm depth below the ground surface and +2.5 cm above the surface.We coated each iButton with silicon sealant to prevent leaking and sensor damage; however, a 25 % sensor failure rate was still experienced.A 10 min sampling interval was used for a minimum of 7 days during July and August of 2013 for each temperature time series.
The installation locations chosen represented the two types of seepage observed with the FO-DTS, and 1-D vertical temperature data were collected synchronously with DTS deployments.Two additional control deployments of 1-D temperature profiles were installed within or below drainage ditches.We assume that under low surficial flow conditions the system is at quasi-steady-state, allowing us to estimate (upward) seepage flux from measured surface water, groundwater, and intermediate-depth temperatures using the analytical solution to the heat transport equation derived by Turcotte and Schubert (1982) and modified by Schmidt et al. (2007).A flux value was calculated for each collected data time step, and was averaged for each profile for the final reported flux value.Flux values were calculated 4 times for each profile using the range of peat porosity and range of thermal conductivity values.The thermal parameters utilized for the 1-D heat transport equation are shown on Table 2.

Assessment of environmental isotopes to infer groundwater flow paths
To trace the source of the groundwater flow paths contributing to discharge, we use δ 18 O and δ 2 H to distinguish between local recharge (short flow paths) and regional recharge (long flow paths).Turcotte and Schubert (1982) and modified by Schmidt et al. (2007).K s is the thermal conductivity of the solid, K f is the thermal conductivity of the fluid, and n is the porosity of the matrix.The density of the fluid and heat capacity of the fluid multiplied together are the volumetric heat capacity of the fluid (ρ f c f , J m −3 K −1 ).water (seasonally), groundwater seepage (August 2013), and pore waters (October 2013).The four pore water samples were acquired through a manual press of samples from Russian peat cores 0-1 m below the ground surface and subsequently filtered for analysis.
Values of δ 2 H-H 2 O and δ 18 O-H 2 O were measured by wavelength scanned cavity ring-down spectrometry on unacidified samples with a Picarro L-1102i WS-CRDS analyzer (Picarro, Sunnyvale, CA).Samples were vaporized at 110 • C. International reference standards (IAEA, Vienna, Austria) were used to calibrate the instrument to the Vienna Standard Mean Ocean Water and Standard Light Antarctic Precipitation (VSMOW-SLAP) scale and working standards were used with each analytical run.Three standards that isotopically bracket the sample values are run alternately with the samples.Secondary laboratory reference waters (from Boulder, Colorado; Tallahassee, Florida; and Amherst, Massachusetts) were calibrated with Greenland Ice Sheet Precipitation (GISP), SLAP, and VSMOW from the IAEA.The isotopic composition results use a rolling calibration, which calculates each sample's error by the three standards run closest in time to the sample.Long-term averages of internal laboratory-standard analytical results yield an instrumental precision of 0.51 ‰ for δ 2 H-H 2 O and 0.08 ‰ for δ 18 O-H 2 O.
The USGS wells were sampled for groundwater isotopic compositions within the PCKD aquifer, providing regional groundwater values for the aquifer and defining the expected annual range of isotopic values for local precipitation (Table 1).The regional groundwater trend line was generated by fitting a linear regression through the USGS well isotope data from the regional PCKD aquifer.

Results
As an initial evaluation of the groundwater contribution to the site, we conducted differential discharge gauging measurements on 15 September 2013.The locations of these measurements are indicated by the purple circles in Fig. 1b.The stream gained 6 L s −1 discharge through Cell 7 from the Arm Pond input to the confluence with Beaver Dam Brook (1.5 km), equal to an average of 0.004 L s −1 per meter of river length (Fig. 1b).Cell 3 and 4 gained 113 L s −1 from the Beaver Dam Pond input to the confluence with the east side river (1 km), equal to an average of 0.113 L s −1 per meter of river length.At other wetland sites seepage flux magnitudes and directions have shown to be temporally transient (Fraser et al., 2001;Sebestyen and Schneider, 2001); however, due to the consistent high hydraulic gradient in the regional aquifer and the small watershed, we assume that temporal dynamics are insignificant within our data set and sufficiently static to describe the present-day conditions.This assumption is supported by the two seasonally distinct infrared surveys resulting with similar seepage distribution results.

Resolving peatland basin structure
The interpolation of the basal surface, or the peat-sand contact beneath the peat from GPR data, indicates four isolated peat depressions at the site, two depressions in Cell 6 and Cell 7 and two in Cell 3 and Cell 4. Cells 6 and 7 have a maximum peat thickness of ∼ 7 m and a gradual curvature of the peat-sand interface than the western cells, Cell 3 and Cell 4 (Fig. 2).The western cells show a maximum peat thickness of ∼ 10 m, and relatively high curvature values.The basin structure of the western cells is also more complex than Cell 6 and 7, as Cell 3 and Cell 4 have pronounced undulations in the basal peat-sand contact surface, creating dramatic changes in basin shape.Particularly, there is a notably high curvature of the basal peat-sand interface along the western edge approximately 30 m from the margin.The GPR profiles illustrate multiple series of normal faults beneath the peat body that are consistent with ice melt-out and/or collapse features (Fig. 2c) typical of kettle-hole origin (Kruger et al., 2009).

Thermal evaluation of groundwater seepage
Surface water temperatures in the main channel and ambient drainage ditch environments generally show high standard deviation, indicative of a coupling between these surface waters and air temperatures, and mean water temperatures closely tied to the seasonal surface temperature average, also indicative of surface water dominance.FO-DTS surveys were designed to detect low standard deviation and consistent mean temperature anomalies from these background conditions, which is indicative of groundwater inflows.The temperature results of both these surveys are presented in Hare et al. (2015).Results from both TIR and FO-DTS identified two categories of thermal anomalies: type 1 anomalies manifest as temperatures with relatively low standard deviation through time, and an anomalous heat signature that is seasonally warmer or cooler than regional groundwater temperature by approximately ±3-5 • C; and type 2 thermal anomalies also have a low standard deviation, but tempera-tures more closely resemble regional groundwater temperatures (10-11 • C). Figure 3 shows time series data collected with the FO-DTS and illustrates each of the major thermal signatures shown on-site: temperatures of groundwater, the main channel, a drainage ditch, and the two thermal anomalies.We interpret these two anomalies to correspond to two modes of seepage, type 1 thermal anomalies correspond to matrix seepage, and type 2 thermal anomalies correspond to preferential flow-path (PFP) seepage.The two seepage types are clearly differentiated through thermal signatures, and can be isolated using the average and standard deviation of temperatures with time.The TIR surveys also revealed these two distinct types of seepage, which were present in both the summer and winter surveys (Fig. 4).
TIR surveys and FO-DTS data indicate that most groundwater input likely occurs along the western edge of the Cell 3 and Cell 4, where peat is thinner or where there is strong sand-peat contact curvature in peat basin shape (Fig. 5).Isolated locations of consistent temperatures similar to groundwater temperatures and anomalously low standard deviations exist along the linear location of highest peat-sand contact curvature near the western edge of the cells, as well as along edge areas with the thinnest peat.The isolated, unique locations of PFP seepage that occur within the deeper peat represent a distinct seepage process from matrix seepage and PFP seeps along the edge of the peat.
During the March infrared survey, a high density of ∼ 1-5 cm diameter flowing macropores within the peat was discovered in Cell 3. The water discharging from these macropores exhibited groundwater seepage temperatures (Fig. 6) and led us to term this mode of PFP seepage.This observation is similar to the peat macropores or "peat pipes" described in previous peatland research (e.g., Briggs et al., 2016;Cunliffe et al., 2013;Holden, 2004;Smart et al., 2012;Vandenbohede et al., 2014), but the concentration of macropores in this singular location makes the northwest cell macropores observation unique.We measured high 3.0 L min −1 flux PFP seeps with a seepage meter.Despite the very few locations of PFPs, their high fluxes have the potential to contribute significantly to the groundwater gain across the site (Poulsen et al., 2015).The peat thickness map (Fig. 5) indicates that the zone of high macropore density is an area of peat thinning reaching a minimum peat thickness of 3 m, and also a location of high curvature (center of cell 3).Rossi et al. (2012) describes similar correlation to peat thinning at a site in Finland.

1-D vertical temperature profiles
The two seepage types and two ambient drainage ditch locations were monitored with 1-D vertical temperature profiles for 7 to 10 days.We expected to observe significant upwelling at this site, which we could easily identify by a rapid attenuation of the diurnal signal with depth coupled with a characteristic convex upward shape of mean temper- Figure 3. Fiber-optic distributed temperature sensing (FO-DTS) temperature time series from four 1 m segments of cable to illustrate the characteristic thermal signatures at Tidmarsh Farms.The greatest amplitude and variability occurs in the drainage ditches with little flow and significant solar heating (red), followed by the main channel of Beaver Dam Brook (green).Two seepage types are also plotted over 2.5 days: matrix (type 1) seepage, with very low variability (low standard deviation) over time and a mean temperature a few degrees higher than groundwater (light blue) and preferential flow-path (type 2) seepage with a mean temperature nearly equal to groundwater (dark blue).
ature with depth (e.g., Schmidt et al., 2007).Temperatures from all four 1-D vertical temperature profiles are distinct from one another; however, all the temperature profiles, including the "ambient" drainage ditches, are consistent with upwelling of groundwater (convex upward shape of mean temperature with depth in Fig. 7).The surface temperature of the ambient drainage ditches (temperature profiles 3 and 4) is similar to the diurnal temperature cycles measured with FO-DTS, and were used as background data for the heat signature of the site.The 1-D fluid flux calculations of the temperature time series of the two drainage ditch locations yielded a range of −0.028 to −0.031 and −0.067 to −0.074 m d −1 .
Temperature profiler 1 was installed at a location with a surficial temperature of 13-14 • C in August 2013.The total peat thickness at this location is 50 cm, and, consistent with groundwater upwelling, minimal diurnal signal propagates to depth, and surface water exhibits relatively low variance in temperatures over time.Thermal time series estimates of flux show a modest −0.146 to −0.163 m d −1 upwelling through the peat at this seepage location.
Finally, temperature profiler 2 was installed in a location with a surficial temperature consistent with groundwater temperatures of 10-11 • C in August 2013, and temperatures with depth exhibit a groundwater thermal signal throughout the entire profile.Even close to the bed interface, the streambed thermistor (2.5 cm) shows slight thermal shifts (σ = 0.096 • C), which are near to the resolution of the instrument (0.0625 • C).This unique temperature profile is indicative of high upward flux rates, as the diurnal signal cannot be resolved and there is essentially no downward conduction from above; therefore, we were unable to use the steadystate analytical solution to estimate a flux rate.However, in Hydrol.Earth Syst. Sci., 21, 6031-6048, 2017 www.hydrol-earth-syst-sci.net/21/6031/2017/ July 2015, we deployed a seepage meter at this location and measured fluxes in excess of 3 m d −1 , rates which exceed the limits for analytical flux calculations.

Groundwater discharge source areas
Groundwater discharge to the wetland complex is a mixture of shallow and deep regional flow paths.Isotopic analyses of waters from wells in the upgradient portion of the PCKD aquifer (blue circles in Fig. 8) fall along a regional groundwater trend line.We interpret this regional trend line to be characteristic of the annual isotopic composition of recharge water to the region as well as local groundwater recharge in the topographic watershed of Tidmarsh.These upgradient groundwater isotopic values plot left of the global meteoric water line (GWML) (Craig, 1961), which reflects local and regional vapor recycling and a characteristic mixture of vapor sources (Koster et al., 1993).The one exception to this line is the USGS well MA-PWW 494 in Plymouth, MA which is similar to Tidmarsh in that it is downgradient of the recharge area of the PCKD aquifer.This water falls to the right of the regional groundwater trend line.Discharging and shallow groundwaters at the wetland site plot close to but off the regional groundwater trend line.The blue diamonds (Fig. 8) represent a monthly sampling of wetland surface waters that depict a significant clustering to the right of the regional groundwater trend and evolve along a line tangential to that intersecting the deep TM groundwater.Uncharacteristically, the deepest sampled groundwater at the site (> 15 m) falls to the right of the GMWL (orange circle), suggesting this water has experienced a significant enrichment in the heavy isotopes due to evaporation processes.Repeated sampling of this water reveals a consistent isotopic composition that suggests the deep groundwater beneath Tidmarsh is isotopically enriched due to evaporation from open-water bodies in upgradient kettle ponds.The headwater seepage area and the strong discharge seepage area (large pink and red triangles in Fig. 8) in the interior of the wetland complex fall along a line that represents either a mixture of this evaporated water and the regional groundwater trend (finely dashed line) or itself is simply evaporatively evolved water.Both interpretations suggest that the source of water to the shallow groundwater wells and the large volume springs in the interior of the wetland complex are distinct.This indicates that the local flow path from the southwest to the northeast is the large-scale hydraulic gradient that dominates the observed seepage patterns.The orientation of peatland basin slope break and the regional groundwater gradient also intercepts the southwest corner of the peatland where numerous high-flux groundwater seeps are located.

Groundwater discharge types
Two types of groundwater discharge (or seepage) were identified using thermal methods, as detailed in Sect.3.2: PFP discharge areas that have regional groundwater temperature (e.g., 10-11 • C) and matrix seepage locations with ground-water discharge at temperatures that are offset (±3-5 • C) from regional groundwater temperature but have very low variance compared to expected diurnal variations and are also significantly distinct from local surface water temperatures.Both seepage types appear to strongly buffer stream temperatures, illustrated by low variance when examined through time (FO-DTS data).A low variance could have also been caused by mobile sediment (Sebok et al., 2015); how- the Tidmarsh Farms area surface water (diamonds), groundwater (circles), shallow well, and seep sources (triangles).The regional groundwater trend line was derived from samples from relatively shallow, regional USGS wells (blue dots), consistent with a relatively humid climate at the site.
ever, within this peatland environment this process is not expected, nor was it observed.The identification of these two distinct seepage types using multiple methods and during distinct seasons indicates different mechanisms for genera-tion of each of these seepage patterns.Figure 5 combines both matrix and PFP seepage observed with either FO-DTS or TIR to evaluate spatial patterning and consistencies, and shows how the two types are related to one another as well as to patterns of high basal curvature.Consistent (low standard deviation) and groundwater-like temperatures (10-11 • C) of the PFP seepage indicate very high flux (> 3 m d −1 was confirmed with seepage meter measurements).Given the low vertical K of peat matrices, sustaining such high fluxes would require seemingly implausible hydraulic gradients, certainly far above the vertical hydraulic gradients observed on-site.Therefore, it is highly likely that this seepage does not occur as flow through the peat matrix, but instead as focused, high-discharge, conduit flow, consistent with "short-circuit discharges" described by Conant Jr. (2004).Focused flow in conduits through the peat was observed in the field at Tidmarsh Farms (Fig. 6), and by Briggs et al. (2016), and has been documented through visual descriptions of peat pipes, or macropores at other locations (Baird, 1997;Beckwith et al., 2003;Cunliffe et al., 2013;Holden, 2004;Smart et al., 2012;Wallage and Holden, 2011).However, the spatial extent of these preferential flow zones has not been previously demonstrated.Due to their high flux, physical isolation, and focused nature, we refer to these types of seepage as PFP seeps.
Data represented by matrix seepage show that surface water diurnal temperatures are also buffered in these zones and are distinct from most ambient surface temperatures.This observation could indicate shallow aquifer groundwater discharge, which is more influenced by atmospheric temperatures than deeper regional flow (Kurylyk et al., 2014a;Menberg et al., 2014).However, consistent temperatures in the site's shallow groundwater wells and 1-D temperature profiles indicate that these seepage temperatures are controlled by a lower flux rather than distinct atmospherically influenced shallow flow paths.These matrix seeps indicate that while vertical upwelling fluxes are present, they are much smaller than PFP discharge zones and must be controlled by a different mechanism.Thermal profilers yielded vertical flux rates consistent with a low to moderate upwelling though porous media according to Conant Jr. ( 2004), which would be typical of the hydraulic properties associated with peat, and thus is the reason we refer to locations with this signature as "matrix" seeps.The two seepage types, PFP and matrix seepage, are similar to the "point" and "diffuse" peat seepage categories defined by Rossi et al. (2012) but, rather than focusing on the area of influence, instead highlight the physical structure that governs the process which ultimately generates seepage in these peatland seepage zones.

Subsurface structural control on the spatial distribution of seepage types
Matrix seeps were plentiful within approximately 30 m of the peatland edge (Fig. 5), consistent with margin seepage observed in lake environments (Rosenberry et al., 2010;Sebestyen and Schneider, 2004;Sebok et al., 2013;Winter, 2001) and other wetlands (Freeze, 1988;Labaugh et al., 1998).The peat is 0.1-3.0m thick along the margin where matrix seepage occurs (Fig. 3), which is generally significantly thinner than locations of observed interior PFP seepage.Matrix seeps generally occur in the thinnest peat zones and typically decrease rapidly with distance from the peatland edge toward the interior slope change, after which no thermally distinct groundwater discharge points are observed (Fig. 5).While evidence for PFP seepage does occur as well in these shallow areas, matrix seepage is more consistent within this shallow peat environment.This is shown as a conceptual model in Fig. 9, based on temperature data collected proximal to GPR line 7.1 (radargram shown in Fig. 2c).Similar landscape-scale observations have been made within lakes and wetlands (e.g., Cherkauer and Zager, 1989;Sebok et al., 2013), and as kettle-hole peatlands typically form from initially open-water bodies, there are logical similarities in basic processes between the two environments.Discrete seepage zones may reflect zones of higher effective K than the surrounding peat matrix, which could be explained by littoral-zone migration in the lake to wetland evolution as the water table fluctuates and migrates.In lake environments, diffuse matrix seepage occurs because of an increase in K at the edge of the lake caused by "erosional deposition," whereby focused wave and current action dis- rupt and erode sediments, particularly mobilizing the finest sediments elsewhere and concentrating larger particles, indicative of these higher-energy environments in these locations.Preferentially stronger flow paths are thus concentrated at the break in land surface slope (Blume et al., 2013;Casson et al., 2010;Cherkauer and McKereghan, 1991;McBride and Pfannkuch, 1975;Rosenberry et al., 2010;Winter, 1981).Previous work proposes that seepage flux decreases exponentially with distance from the shore of a lake (Cherkauer and Zager, 1989;McBride and Pfannkuch, 1975) (Newby et al., 2000(Newby et al., , 2009)).Therefore, we hypothesize that the extent of the matrix seepage observed along the western edge of the peatland is a result of this lake transgression and coincident decrease in deposition of organic material.Here the lower K of the peat matrix intersects with shallow groundwater flow paths, strongly affecting lateral hydraulic gradients and driving upward flux; a process which likely generates much of the observed matrix seepage (Fig. 9).This observation is supported by similar seepage processes observed in riverine systems (Sophocleous, 2002), wetland (Larsen et al., 2007), lake (Bakker and Anderson, 2002;Winter, 1981), and hillslope environments (Shaw et al., 2017;Winter et al., 1998).
In contrast to the matrix seepage, PFP seepage was less common and spatially disconnected from similar flux seeps (Fig. 5).Similar to matrix seepage, PFP seepage exhibits low standard deviation of temperature (Fig. 3), but PFP seep temperatures were much closer to average regional groundwater temperature.This indicates that PFP seepage waters have very short residence times within peatland sediments, which may have important implications for nutrient transformations within them.At some PFP seeps the peat is generally thicker and located more toward the interior of the peatland rather than along the margin where matrix seepage zones are found in addition to being found between the peatland edge and the area with high basal curvature values (Figs. 5 and 9).Typical interior PFP flow-path lengths from the sandy aquifer below the peat to the surface should be much greater than for matrix seeps; however, the thermal signature seems to contradict this; therefore, PFP seepage zones must be generated through a unique hydraulic process from matrix seeps.Since PFP seeps at Tidmarsh Farms correlate with significant slope changes, or locations of high curvature, these isolated seepage zones must be generated by an abrupt change in horizontal K, and the PFP seep locations closer to the edge may be a result of zones of inherent matrix weaknesses, such as varying degrees of humification caused by vegetative difference and water level, or other disruptions in the peat matrix, including plant rooting and desiccation "cracks" as proposed by Smart et al. (2012) (Fig. 9).
An abrupt change from high to low K has long been known to promote the transition from horizontal to vertical flow (Freeze and Witherspoon, 1967).Lowry et al. (2009) hypothesized this process to explain developed seepage within the interior of a peatland through using 3-D numerical groundwater flow models.As horizontally flowing regional groundwater encounters a low-conductivity peatland, it is forced to go through or around it, causing pressure to increase where the abrupt change in the K from the sand to catotelm peat matrix occurs (Fig. 9).PFP seeps develop as fast pathways to the surface and as pressure-relief valves, where these localized increases in aquifer pressure at the base of the peat matrix translate into strong, sustained discharge of unaltered regional groundwater to the surface.Rosenberry et al. (2010) note that in lake bottoms, a significant upward seepage velocity can maintain a locally high K, as the upward force may suspend smaller particles within the water column.Particulate organic matter and lacustrine sediment have a very low settling velocity, and therefore if the upward force that groundwater seepage induces is greater than the settling velocity, only organic matter with a high mass will be able to accumulate over these lake seepage locations.This would cause the peat matrix to have a relatively high porosity and a high permeability compared to its surrounding very low permeability matrix.These locations will continue to be zones of weakness through the formation of the peatland.Thus, we hypothesize that high-flux PFP seepage zones persist through the transition from lake to peatland environment due to the inability of fine sediments and organic matter to accumulate over these high-flux locations.Still, these locations of consistently high hydraulic gradient will also continually take advantage of inherent matrix weaknesses.However, the underlying mechanics of PFP seepage in the deeper interior peat are caused by the interception of the regional groundwater gradient and high-curvature peat subsurface structure (Fig. 9).
The orientation of peatland basin slope break (high basin profile curvature) and the southwest-to-northeast regional groundwater gradient dictates the observed pattern of strong seepage along the western boundary, which is supported by isotopic analysis.PFP and matrix seep waters both exhibit isotopic signatures consistent with a mixture of local groundwater and regional recharge signature (Fig. 8).This observation is further reinforced by the increase in net groundwater gain through the western cells, as well as a large number of PFP seeps in the southwestern portion of the site (Fig. 5).

Conclusions
Subsurface basin shape exhibits significant control on the spatial distribution of groundwater discharge within peatland environments.As horizontal groundwater flow intercepts the peat matrix, two types of seepage develop: matrix and preferential flow-path seepage.Matrix seepage is defined by a low standard deviation in temperature and surface temperature similar to groundwater ±3-5 • C, consistent with relatively low-flux seepage.Low fluxes are produced where the regional groundwater flow paths intercept the low-K peat at the basin 'shoreline', inducing upward flow through relatively thin (0.1-3.0 m) peat.The second type of observed discharge, PFP seepage, has a surface temperature essentially indistinguishable from deep regional groundwater temperature.This indicates very strong upwelling fluxes at these locations and little time for conductive heat losses or gains.Locations of PFP seeps appear along the periphery of the peatland, but more notably also correlate with high rates of basal peat slope change (curvature) of the peat basin (Fig. 9).These seeps develop where the regional groundwater flow path intercepts a secondary slope change and where there is a stark change in K between the high-K sand aquifer material and the low-K peat.Together, these physical features generate large pressures, induce localized zones of high vertical hydraulic gradient, and drive large seepage fluxes upward.Because PFP seeps typically occur in locations with thicker peat and yet maintain close to groundwater temperatures, they must have a much higher vertical hydraulic gradient and/or higher effective K than the matrix seeps.Through multiple lines of evidence, we conclude that the development and spatial distribution of minerotrophic peatland seepage is strongly controlled by the interactions among the subsurface basin structure, physical processes within the peat structure, and hydraulic gradient.
Through our results, we establish a predictable pattern of seepage, consistent across the coastal site that is explained by knowledge of basin shape and regional hydraulic gradient.This information provides valuable insights for water resource managers to better understand the natural forces driving groundwater seepage.This knowledge, in turn, may be used in the restoration design of degraded peatland systems.Knowing where seepage is expected to occur naturally across a site allows for the development of more sustainable restoration designs that work with the land, and not against it.In retired cranberry farms, for example, channels may be relocated to intercept springs to maintain cooler water temperatures.This knowledge can also guide the location of targeted intensive grading.For example, as was done at Tidmarsh Farms, the dense cranberry mat can be broken up mechanically to encourage groundwater expression on former dry farm surfaces and access native seed banks below.Incorporating this data into a restoration design will greatly aid the ability to predict and achieve desired ecosystem outcomes, making restoration projects more efficient, both ecologically and monetarily.
This research provides a process-based investigation of the subsurface hydrodynamics within a peatland.While a peat matrix exhibits strongly heterogeneous and anisotropic tendencies, large-scale patterns occur and can be predicted.These patterns are dependent on basin shape, peat accumulation history, and underlying aquifer flow paths.The importance of groundwater flow paths surrounding the peatland and resulting seepage patterns emphasizes that peatlands are not isolated entities from the groundwater system and cannot be treated as such.Observed large-scale seepage patterning provides insight that may help explain vegetation patterning, macropore development, and other localized peat dynamics that have been unidentified in the past, and greatly aid peatland management and restoration to establish more naturally sustainable, efficient practices.
Data availability.The underlying research data from this manuscript can be accessed through email request to the corresponding author.
Competing interests.The authors declare that they have no conflict of interest.
Disclaimer.The views and opinions expressed in this article are those of the authors and do not necessarily reflect the official policy or position of AECOM Technical Services.
illustrates the location of field measurements.Dif-

Figure 2 .
Figure 2. (a) Map of total peat thickness beneath Tidmarsh Farms based on GPR data.GPR data collected along linear transects shown here (black lines; pink lines in Fig. 1b) were interpolated and contoured to show peat thickness (colors) on the 2-D surface map.Zones of medium and high curvature (the second derivative of the thickness) of the peat-sand interface are shown as grey and black pixels, respectively.(c,d) Three example cross sectional profiles, or radargrams, illustrate a distinct reflector at the basal peat-sand contact.Peat is shaded red.Sediment cores shown as yellow lines (hand cores) and orange lines (vibracore) were used to constrain the GPR velocity data.High curvature is highlighted in green boxes.

Figure 4 .
Figure 4. Thermal infrared (TIR) images recorded on 30-31 July 2013 (Summer) and 21 March 2014(Winter)  at Tidmarsh Farms.Visiblelight images are shown in the bottom left of March images, but not July, as these surveys were conducted at night to limit issues associated with reflectance.TIR images illustrate the two types of seepage in both seasons: type 1 preferential flow-path seepage that is characterized by discrete discharge points very close to groundwater temperature with high flux, and type 2 matrix seeps that are diffuse and 3-5 • C warmer or cooler than groundwater and lower flux.

Figure 5 .
Figure 5. Map of seepage at Tidmarsh Farms determined with fiber-optic distributed temperature sensing (FO-DTS, squares) and thermal infrared (TIR) surveys (circles).Background shaded region(s) match the bounded area from Fig. 1b, and darker background shading delineates zones of high curvature (the 2nd derivative of the thickness) of the peat-sand interface (Fig. 2).For both methods, light purple to pink symbols indicate matrix (type 1) seepage, and dark blue indicates locations of PFP (type 2) seepage.From FO-DTS data, a location was tagged as seepage if the standard deviation was less than 1.5 and the temperature was less than 15 • C for the matrix and less than or equal to 11 • C for PFP seepage.From TIR surveys, seepage was distinguished by temperatures of 9-11 • C for interior seepage, and 11-15 • C for matrix seepage.The location of GPR line 7.1 is shown on this figure to reference data for the conceptual model in Fig. 9.

Figure 6 .
Figure 6.Thermal infrared (TIR) image from 21 March 2014 at Tidmarsh Farms illustrating PFP (type 2) seepage.Many macropores are observed in both the infrared (slightly smaller) and the visual image.These seeps are located in the middle of cell 3 (Fig. 1b), where peat is ∼ 3 m thick and dramatically thinning.

Figure 7 .Figure 8 .
Figure7.Temperature profiles vs. depth at Tidmarsh Farms recorded in 30-31 July 2013.For each profile, the range of air temperatures and groundwater (GW) temperatures are shown as bands of pink (air) and dark blue (GW).Locations 1 and 2 (profiles 1 and 2) show the influence of upwelling GW, expressed as type 1 preferential flow-path (PFP) seepage (profile 2) and type 2 matrix seepage (profile 1).The convex upward shape of temperature-depth profiles 3 and 4 is also consistent with upwelling GW.

Figure 9 .
Figure 9. Conceptual model illustrating the mechanism for development of matrix (type 1) seepage (pink arrows) and preferential flow path (PFP) or interior (type 2) seepage (blue arrows), corresponding to locations in winter TIR images.Thick black lines represent groundwater flow direction, and the yellow-green box indicates the region of high basin curvature.The brown basin represents peat in a typical basin shape based on GPR line 7.1 (Fig.2c).Conceptual matrix and PFP seepage locations are based on the temperature data recorded proximal to GPR line 7.1 and winter TIR images from this same transect.PFP seeps found in the thicker peat are associated with locations of high basin curvature where strong vertical gradients drive focused, higher-flux seepage through pre-existing weaknesses in the peat matrix.

Table 1 .
USGS groundwater wells δ2H-H 2 O and δ18O-H 2 O isotopic data used to establish the regional groundwater trend.