Aquifer-Peatland Hydrological Connectivity and Controlling Factors in Boreal Peatlands

The conditions in which groundwater inflow occurs in boreal peatlands and its contribution to peatland water balance are still poorly understood. The objectives of this research were to quantify the hydrological connectivity between a surficial aquifer and a peatland, and to identify the controlling factors in boreal peatlands of north-central Quebec (Canada). The peatlands were instrumented with piezometers and groundwater levels were monitored during two growing seasons. Hydraulic conductivities were measured on peat cores and in situ, groundwater inflows and outflows were calculated using the Darcy equation. The peatland water budgets were simulated for the two peatlands with a steady-state groundwater flow model to verify flow hypotheses, to quantify unmeasured flows and to explore recharge scenarios leading to changes in connectivity. The two peatlands have contrasted water budgets, with recharge representing the largest inflow (78%) and subsurface runoff representing the largest outflow (85%) the peatland with the smallest catchment area (Misask). The peatland with the largest catchment area (Cheinu) is also located downgradient within the regional watershed. Its inflows are dominated by groundwater (56%) and its outflows are mostly towards subsurface runoff (74%). The two peatlands are in conditions of precipitation excess and a recharge reduction would not affect their peatland heads markedly (<10 cm). However, recharge changes could induce larger modifications in groundwater inflows and outflows for the peatland with a larger catchment area. The dominating peatland hydrological functions are thus contrasted at the two sites, and it is hypothesized that the water table depths thresholds triggering changes between storage, transmission and runoff functions are also different. Although further studies remain to be done to understand how hydrological conditions change through time, and ultimately what are the long-term impacts of a changing climate on hydrology, vegetation and carbon accumulation, this work shows that understanding peatland hydrology requires to consider hydrological conditions beyond the peatland limits.


INTRODUCTION
Peatlands cover 12-17% of Canada (Tarnocai et al., 2000) and approximately 9-12% of the province of Quebec (Payette and Rochefort, 2001). Canadian boreal peatlands correspond to almost half of the Earth's wetland cover (Richardson et al., 2021). In the north-central region of Quebec, peatlands are mostly characterized by poor fens dominated by a succession of elongated strings and pool microforms parallel to the slope (Payette and Rochefort, 2001;Robitaille et al., 2021) comparable to Fennoscandia and northwestern Russia boreal aapa mires (Seppälä and Koutaniemi, 1985;Foster and Fritz, 1987).
Peatland hydrology has received much scientific attention over the past decades (Whittington and Price, 2006;Spence et al., 2011;Waddington et al., 2015;Goodbrand et al., 2019). Their carbon storage capacity is another focus of many studies (Yu, 2011;Gallego-Sala et al., 2018;Loisel et al., 2021). Peatland resilience to natural and anthropic disturbances such as drainage (e.g., Menberu et al., 2018;Harris et al., 2020), roads (Saraswati et al., 2020;Elmes et al., 2021), fires (Bourgeau-Chavez et al., 2020;Nelson et al., 2021), and climate change is also of growing concern. According to some climate projections, northern peatlands will experience an increase in vapor pressure deficit (Yuan et al., 2019), resulting in increased evapotranspiration that could exceed that of the surrounding boreal forest (Helbig et al., 2020). However in Eastern Canada, and more specifically in northern Quebec (Canada), climate projections for the middle of the 21st century predict a warmer and wetter climate (Ouranos, 2015). Generally speaking, these conditions can lead to important changes in groundwater recharge (Ireson et al., 2015;Schneider et al., 2016), which can affect water table depths, and consequently impact the influx of groundwater to peatlands (e.g., Hokanson et al., 2020).
Although peatland hydrological functions have been the subject of few quantitative studies, it has been shown that peatlands can recharge the aquifer, store and transmit water, and sustain discharge to rivers and streams (cf. Goodbrand et al., 2019). It is recognized that these flows can be highly variable depending on the geological and climatic contexts . Groundwater flow models have been used to quantify flows exchanged between aquifers and peatlands, showing that groundwater inflow can represent up to 50% of the peatland water budget in southern Quebec (Bourgault et al., 2014;Levison et al., 2014). However, aquifer-peatland hydrological connectivity is still poorly documented, especially in boreal peatlands that have been relatively seldomly studied in northern Canada. Sensitivity of aquifer-peatland connectivity to recharge has been addressed in the scientific literature with groundwater flow models that explicitly represent the hydrogeological characteristics of both the geological media and the organic deposits (e.g., Bourgault et al., 2014;Levison et al., 2014;Rossi et al., 2014;Thompson et al., 2015;Li et al., 2019). Autio et al. (2020) summarize studies that have used integrated surface-subsurface flow models including peatlands to study this connectivity. These authors report that peat properties and landscape topography are crucial parameters in models aimed at quantifying interactions between surface and subsurface reservoirs.
The objectives of this research were to quantify aquiferpeatland hydrological connectivity and identify controlling factors in boreal peatlands of north-central Quebec (Canada). Groundwater levels were monitored over two growing seasons, peat and aquifer hydraulic conductivities were measured in situ and in laboratory conditions, groundwater inflows and outflows were calculated using field measured data. Groundwater flow models were built to verify flow hypotheses, to quantify unmeasured flows and to explore recharge scenarios leading to changes in connectivity.

Geology
The two studied peatlands are located in the north-central region of Quebec (eastern Canada), near the Stornoway Inc. diamond mine (named Renard mine). The retreat of the Laurentide Ice Sheet occurred between 7,800 and 7,450 cal BP (Dyke et al., 2003). Surface deposits are dominated by sandy to silty (including boulders) glacial till showing elongated erosional (drumlins) forms associated with subparallel eskers oriented at N220° (Lamarche and Hébert, 2020). Available drilling data show that glacial till can reach 30 m (Desbiens, 2008) and glaciofluvial deposits can exceed 40 m (Lamarche and Hébert, 2020). Surface deposits are highly permeable and there are no marine or lacustrine deposits. The vertical stratigraphy is known precisely, but considering the post-glacial environment, it is probably not highly complex and does not lead to large vertical heterogeneity. The hydraulic properties of Quaternary deposits and underlying bedrock have not been measured prior to this study. The surficial aquifer is unconfined, and the water table is expected to be shallow throughout the region, as indicated by data from a 5 m-deep well in the granular aquifer (aquifer depth unknown; average water table depth 35 cm).

Misask and Cheinu Peatlands
Peatlands in the region developed in depressions of the Precambrian Shield. The Misask (52°43′ 27″ N, 72°12′ 50″ W) and Cheinu (52°38′ 48″ N, 72°11′ 31″ W) peatlands (unofficial names) are located at the ecotone between the open (taiga) and closed boreal forest biomes ( Figure 1). Peat accumulation initiated ca. 6,430 and 6,560 cal BP (Robitaille et al., 2021), suggesting early post-glacial initial vegetation cover. This region represents the biogeographic limit of ombrotrophic peatlands, which dominate the boreal region and are gradually replaced by oligotrophic patterned fens with abundant pools moving northwards (Payette and Rochefort, 2001). Both peatlands show a minerotrophic surface pattern with alternating elongated strings and pools parallel to the slope, and have undergone recent ombrotrophication over the past decades (Robitaille et al., 2021). Surface vegetation is generally dominated by Sphagnum moss species, with local changes controlled by microtopography with herbaceous species (Carex spp. and Trichophorum cespitosum) mainly found in depressions while ericaceous species (Chamaedaphne calyculata and Kalmia angustifolia) are found in strings (Robitaille et al., 2021). Peripheral minerotrophic areas are characterized by thinner peat (<1 m), with black spruce (Picea mariana) and larch (Larix laricina) and show recent and ongoing paludification.
A detailed digital elevation model (DEM) was produced using photogrammetric data from a drone survey (Stornoway Mine, unpublished data). Elevations within the peatlands and for all the piezometers were measured using a differential global positioning system (DGPS). The vertical precision of the DGPS data is estimated to be 1 cm in static mode (USGS, 2021). The photogrammetric DTM was geo-referenced to match the DGPS survey with respective XY error and Z error of 0.004 and 0.005 m (Misask) and 0.005 and 0.006 (Cheinu).
Data show that the studied peatlands developed in local depressions created by the drumlin and moraine landscape. These depressions are underlain by boulders and sand-filled interstices inherited from excessive till leaching (Dionne, 1978), which represent an unusually permeable substratum for peatland development. The Misask peatland (0.102 km 2 ) developed in a relatively flat landscape (immediate uphill slope 0.5%) and the peat surface has an average slope of 1%. Shallow pools covering approximately 10% of its total area, are found exclusively in the eastern portion of the peatland. The Cheinu peatland (0.147 km 2 ) developed in a deeper bowl-shaped depression located at the base of a 6.7% slope. The peatland surface has an average slope of 0.5%. The Cheinu peatland shows well-developed pools, perpendicular to the slope and covering 60% of the area. Robitaille et al. (2021) reported a mean peat thickness of 0.71 and 0.87 m for Misask and Cheinu respectively. Maximum peat thickness is reached in the central patterned area of the Misask site (2.08 m) and in the southwest floating mossmat area of the Cheinu site (3.28 m) ( Table 1).
Located in the upper part of the large Eastmain River watershed (46,000 km 2 ), the peatlands are also included within the Misask River watershed, an important tributary of the  Eastmain River. Their respective catchment areas cover 0.226 km 2 (Misask) and 0.487 km 2 (Cheinu), ranging in altitude between 483 and 468 m for Misask (average slope 1.9%) and 482 and 455 m for Cheinu (average slope 3.1%). The peatland watersheds show successions of small kettle lakes connected to minor streams. Permanent outflows are present at the two peatlands in the shallower paludified areas dominated by bryophytes and black spruce. The eastern outlet of the Misask site feeds a small intermittent stream and a second outlet also runs south ( Figure 1). Narrow paludified areas in the southern portion of the Cheinu peatland contribute to a 1.8 m-wide outlet stream that flows into the downstream lake.

Meteorological Conditions
Meteorological data were obtained from a weather station located approximately 15 km north, at the Stornoway mine. Between 2017 and 2019, annual precipitation (P) varied between 779 and 821 mm yr −1 (mean 797 mm yr −1 ). The average annual temperature was −3.0°C, with snow occurring between mid-October and May. The potential evapotranspiration rate (PET) was calculated using the equation of Oudin et al. (2005), based on the mean daily air temperature, and the estimate of Morton (1983) of extraterrestrial radiation. Annual PET for 2017 to 2019 ranged between 341 and 372 mm yr −1 (mean 360 mm yr −1 ). Documented evapotranspiration rates in northern peatlands have been found to be equivalent to or slightly lower than PET (e.g., Wu et al., 2010).

Peatland Instrumentation and Monitoring
Between June 2018 and October 2019, the peatlands were instrumented with five piezometers set in the mineral deposit (P m or P c , "m" stands for Misask and "c" stands for Cheinu; Figure 1) and five observation wells within the peatlands (W m or W c ; Figure 1), each equipped with water level loggers (Solinst) recording hourly. The piezometers consisted of 25 mm-diameter PVC tubes with 0.3 m-long screens at their base. The wells were made of similar pipes but were screened from top to bottom within the peat deposits. A barometer (Solinst) was installed at the Cheinu site for barometric compensation. In the Misask and Cheinu peatlands respectively, five and eight additional noninstrumented piezometers and wells were used for manual head measurements during each field visit (June, July, and September 2018, and June and September 2019). At each site, a metal rod was inserted through the peat (lawn) and fixed into the underlying mineral deposit to estimate peat contraction and expansion between field campaigns. The variation in measured seasonal peat elevations was negligible, with a maximum of 1 cm (Misask) and 2.5 cm (Cheinu).

Hydraulic Conductivities
Two marginal and one central 1 m peat cores were collected at the deepest locations of each peatland ( Figure 1) using a Box corer (Jeglum et al., 1992). The cores were wrapped and brought to the laboratory, where they were kept at 4°C until further analysis. The cores were later cut with a serrated knife into 8 cm × 8 cm x 8 cm cubes. Vertical and horizontal peat hydraulic conductivities were measured using the Modified Cube Method (MCM; Beckwith et al., 2003), following the instructions of Rosa and Larocque (2008).
Slug tests were performed in the piezometers located in the surrounding glacial till and in the fluvioglacial deposits to estimate their hydraulic conductivities. Falling-head tests were conducted by adding a volume of water (i.e., from 10 to 90 ml, depending on recovery time), or by using the logger as a slug, in the PVC tube with recording frequency set to a 2 s interval. Two or three slug tests were performed for each piezometer and the results were interpreted using the Hvorslev (1951) method.

Peatland Inflows and Outflows
Elevations of all piezometers and wells, of the surrounding water bodies, as well as data from three existing observation wells at the Renard Mine airport, were extracted from the DEM-derived topographical data. Horizontal hydraulic gradients were calculated between neighboring peatland well and mineral piezometer pairs to estimate groundwater flows between the aquifer and the peatland using Darcy's Law and the segment method (Rosenberry and LaBaugh, 2008): where q A-P is the unitary flow exchanged horizontally between the peat and the mineral deposits (mm yr −1 ), K h is the geometric mean horizontal hydraulic conductivity (m s −1 ) between mineral deposits (average value derived from the slug tests) and peat (average value for the top 30 cm of peat derived from MCM peat core measurement), l j is the length of each segment around the peatland (m; colored peat contour sections on Figure 1), d p is either the average peat thickness or the maximum peat thickness (m) (the two values are used to set an interval of possible fluxes), i h is the horizontal hydraulic gradient (m m −1 ) estimated from the average heads available for the piezometer-well pair, and A p is the peatland area (m 2 ). The sign of the hydraulic gradient determines whether there is an inflow from the aquifer to the peatland or an outflow from the peatland to the aquifer. Similarly, vertical hydraulic gradients were calculated between peatland wells and piezometers located within the mineral deposits in the well vicinity: where q v is the unitary flow exchanged vertically between the peat and the mineral deposits (mm yr −1 ), K v is the geometric mean vertical hydraulic conductivity (m s −1 ) between mineral deposits (average value derived from the slug tests) and peat (average value at 92 cm derived from MCM peat core measurement), and i h is the horizontal hydraulic gradient (m m −1 ) estimated from the average heads available at the piezometer-well pair for each segment (cf. colored peat contour sections on Figure 1). Groundwater recharge represents the vertical water inflow to the peatland through its unsaturated zone. It was estimated using annual net precipitation values (recharge = P net = P-PET) between 2017 and 2019, based on the hypothesis that there is limited runoff due to the high hydraulic conductivity of the surface deposits. P net varies between 414 and 450 mm yr −1 (average 434 mm yr −1 ). The surface outflow from the peatland at the main outlet of the Cheinu peatland was measured manually during each site visit using a current flow meter (Swoffer). Flow rates were calculated using the velocity-area method (Herschy, 1993). The outlets of the Misask peatland were too diffuse to allow similar field measurement.

Groundwater Flow Model
A finite-difference steady-state groundwater flow model was developed in MODFLOW (Harbaugh, 2005) for each peatland. The MODFLOW platform in steady-state conditions was justified by the available data describing the simulated system (geometry and some hydraulic conductivities values). A grid with a maximum cell size of 128 m was gradually refined to 4 m cells near and within the peatland. Fourteen vertical layers reaching a total depth of 70 m were used, with thickness gradually increasing Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 835817 5 from 0.25 m in the top 2-30 m. The impact of cell size on simulations were tested prior to reaching these values (not shown) and are similar those reported by Levison et al. (2014).
The Misask and Cheinu models cover 15.37 and 7.14 km 2 respectively (approximately 47 and 11 times their watershed areas) to ensure that the boundary conditions have a limited effect on modeling results. Model boundaries were composed of rivers and lakes as full-depth constant head boundaries, connected by no-flow segments identified from the piezometric maps. The lower limit is set as a no-flow boundary. Due to the lack of information on their connectivity to the surficial aquifer, the small lakes located inside the models were represented with constant heads constrained to the top layers ( Figure 2). The impact of the boundary conditions on simulated heads and flows was tested and estimated to be negligible (results not shown). The peat thickness was defined based on data presented in Robitaille et al. (2021).
Outside the peatlands, infiltration zones and hydraulic conductivity zones were established based on the Quaternary deposits map. Each model was discretized into five distinct zones including the peatland ( Figure 2). Sediment types were different for the two peatlands: the Misask watershed includes a more rugged terrain covered with thin glacial till (<1 m) while the Cheinu model includes significant boulder fields with glacial till in its northern portion and fluvioglacial deposits in its southern portion.
Within the peat deposits, hydraulic conductivities for the first meter (layers 1-9) were set to the mean MCM-measured K h and K v values. Hydraulic conductivities for the mineral deposits (glacial till and fluvioglacial deposits) were estimated from slug test results and literature values. A vertical anisotropy (K h /K v ) of 10 was imposed on the peat layer below 1 m and on the mineral deposits.
The Unsaturated Zone Flow Package (UZF) (Niswonger et al., 2006) was used to redirect infiltration to groundwater recharge and seepage. Seepage corresponds to groundwater that outflows to the surface when the water table reaches the surface and is exported out of the peatland. It offers a means to represent groundwater seepage relatively easily within a peatland. The UZF parameter SURFDEP (Niswonger et al., 2006), which represents the maximum amplitude of the topographic variation within a grid cell, was fixed at 0.8 m following the work of Feinstein et al. (2020) who showed that a fen peatland model was not sensitive to this parameter.
Mineral K h and recharge rates were calibrated using PEST (Doherty, 2005) to reproduce average measured heads in the peatlands and in the aquifer. Parameters were then adjusted manually to ensure decreasing values of K with depth.

Recharge Scenarios
Because groundwater recharge is a highly imprecise variable in any aquifer study, the response of heads and flow rates (peatland inflows and outflows) to changes in recharge, as indicators of connectivity response, was performed for the two peatlands. As the peatlands are set in a pristine environment, the main reason for a change in recharge is modifications in precipitation, temperature, and evapotranspiration as a result of climate change. Few studies have estimated groundwater recharge in northern Canada in a climate change context. However, recent studies for more southern temperate locations suggest a range of possible climate-driven changes in recharge, from decreases reaching −50% (Bourgault et al., 2014;Wang et al., 2014), to increases reaching +53% (Scibek and Allen, 2006;Jyrkama and Sykes, 2007). A variety of changes, from −59% to +15%, have been reported elsewhere (Croteau et al., 2010;Levison et al., 2014;Dubois et al., 2022). Here, model reaction to climate-driven changes in recharge was tested with scenarios of −50%, −20% and +20%. The emphasis was put on recharge reduction based on observations from Robitaille et al. (2021Robitaille et al. ( ) who identified (1950Robitaille et al. ( -2017 an increase in temperature and in growing degree FIGURE 3 | Minimum, maximum, and average measured hydraulic conductivities of the three cores for (A) the Misask peatland and (B) the Cheinu peatland. The left panels illustrate vertical hydraulic conductivities (K v ) and horizontal hydraulic conductivities (K h ) for the peat (estimated with the MCM method) and for mineral deposits (estimated with slug tests). Panels on the right illustrate K h /K v ratios.
Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 835817 days, combined with a decrease in precipitation. Similar ranges of recharge changes have been used by Bourgault et al. (2014) and Levison et al. (2014).

Hydraulic Conductivities
Slug tests in the glacial till at Misask ( Figure 3A) show K values ranging between 1.2 × 10 -5 m s −1 and 3.3 × 10 -7 m s −1 (average 4.1 × 10 -6 m s −1 ). Slug tests in the fluvioglacial deposits at Cheinu ( Figure 3B) show slightly higher K values, ranging between 1.2 × 10 -4 m s −1 and 7.6 × 10 -6 m s −1 (average 3.8 × 10 -5 m s −1 ). Because there were no slug tests conducted in washed out till and thin till, their K values are assumed to be equal to those estimated for fluvioglacial deposits and glacial till deposits respectively. At both sites, MCM results showed generally decreasing trends in average K h and K v with depth, with maximum K values observed between the surface and approximately 40 cm depth ( Figure 3). Deeper, K values are relatively constant with depth. The average maximum K h and K v values (topmost cubes) were 8.7 × 10 -4 and 8.6 × 10 -4 m s −1 at Misask ( Figure 3A) and 5.3 × 10 -4 and 5.5 × 10 -4 m s −1 at Cheinu ( Figure 3B) respectively. The top cubes (0-8 cm) of the Cheinu peat cores were dismissed since MCM tests were unable to be properly executed (wax filling peat pores and friability of the living moss). The vertical anisotropy (K h /K v ) ranged from 1.02 (0.04 m) to 7.72 (0.44 m) at the Misask peatland, with low values under 0.20 m followed by irregularly decreasing values deeper in the peat profile ( Figure 3A). A similar trend was observed at the Cheinu peatland, with a top K h /K v of 1.69 followed by a maximum of 10 at 0.20 m followed by a constant decrease with depth and a minimum value of 0.24 at 0.84 m ( Figure 3B).

Water Table Depths and Heads
Monitored wells within the organic deposits at Misask show water table depth (WTD) (Figure 4) variations of between 57 and 6 cm below the surface. Minimum WTD occurred upgradient of the central peatland area immediately following snowmelt (June 2018: W m 05), while maximum WTD was observed in the paludified downgradient area during the growing season (July 2018: W m 08). The Cheinu peatland WTD showed a similar range, from 55 to 7 cm below the surface. Minimum WTD occurred within the central peatland area at the end of the growing season (October 2018: W c 28) and maximum WTD was observed near the woody main outlet during the growing season (August 2018: W c 22). The minimum recorded WTD of 6 cm (Misask) and 7 cm (Cheinu) show that no flooding occurred in the monitored areas between June 2018 to October 2019. At both sites, the most central peatland wells (W m 04 and W c 28) were characterized by shallow water tables, with respective average WTD of 18 and 15 cm.
Except for the central peatland wells (W m 04, shallow WTD) and one peatland well located close to the northern boundary (W m 16, deeper than expected WTD), the Misask peatland generally shows increasing WTD from its northeastern section toward the southwestern area. A similar pattern is visible at the Cheinu peatland. With some exceptions (W c 23 and W c 28), the average WTD generally decreases from north (19 cm) to south (41 cm). These are indications of northeast-southeast and northsouth flow directions within the Misask and Cheinu peatlands respectively. Exceptions to these trends correspond to the forest border where peat is shallow, poorly decomposed with high porosity and associated with current lateral expansion.
In the surface deposits, WTD ( Figure 4) varied from 78 to 6 cm at Misask and from 68 to 1 cm at Cheinu. In both sites, minimum WTD was measured at upstream piezometers following snowmelt in 2018 and the end of the growing season in September 2018. Maximum WTD occurred in downstream piezometers during the growing season, in August 2019 and in July 2019. Piezometers in mineral deposits surrounding the Misask peatland showed generally deeper average WTD compared to those installed at Cheinu. Piezometers at Misask have higher average WTD northwest of the peatland than southwest of the peatland while, the piezometers north of the Cheinu peatland have higher average WTD than those southward.
At the two sites, heads reflect the general groundwater flow directions within the mineral deposits ( Figure 5). At the Misask peatland, heads in the mineral deposits indicate groundwater flows from P m 07 to P m 09, from northeast to southwest and in the direction of the diffuse southwest outlet, and from P m 07 to P m 06, i.e., towards the southern peatland outlet. At the Cheinu peatland, heads in the mineral deposits indicate flow from P c 13 to P c 26 and P c 30, i.e., from north to south. Heads follows a similar direction within the peat, from W c 17 to W c 28 and W c 22, towards the southern peatland outlet.

Hydraulic Gradients and Flows
Hydraulic gradients were calculated for all the sections used in the Darcy flow calculation. A negative horizontal hydraulic gradient inflow to the ecosystem. At both sites, inflows driven by local topography were also found to occur in areas located close to downgradient peatland flow locations. HGGs did not change markedly in time and inversions lasted a maximum of four consecutive days in August and October. At Misask, the largest negative HHGs (absolute value) were observed in the southeastern and eastern portions of the peatland, close to the main diffuse outlet. Results from P m 09-W m 04 (south) and P m 12-W m 04 (east) also indicated that the peatland feeds the aquifer laterally in these areas. At Cheinu, a NHG was only found in the vicinity of the outlet and was markedly smaller than those observed at Misask. Considering that the peatlands developed in shallow depressions, it is interesting to note that topography is clearly not the only factor determining groundwater inflows to the Misask peatland. If that were the case, a positive HHG would be observed along its entire periphery. Topography is apparently a stronger driver of horizontal hydraulic gradients at Cheinu, with only one negative HHG of limited amplitude. Three vertical hydraulic gradients (VHG) were calculated at the Misask peatland. The P m 03-W m 16 pair probably provides the most reliable information, located only 6 m apart, whereas the other two pairs are located 25 m apart. The results showed VHG values of −0.17, −0.10, and 1.16 m m −1 indicating flow from the peatland to the underlying permeable mineral deposits in the northwestern and central-western parts of the ecosystem, and flow from the underling sediments to the peatland in its northeastern part. As for HHG, VHG values were relatively constant over time, with inversions observed over very short periods of only a few days.
Horizontal Darcy flows were calculated for each section using average heads for all the well-piezometer pairs ( Table 2). Darcy inflows of 92 and 269 mm yr −1 were calculated at Misask, and of 111 and 418 mm yr −1 at Cheinu, with the minimum value estimated using average peat thickness and the maximum value estimated using maximum peat thickness at each site. Darcy outflows of 73 and 212 mm yr −1 were calculated at Misask, and of 48 and 131 mm yr −1 at Cheinu, again using average and maximum peat thickness. At Misask, vertical outflows of 30 mm yr −1 (P m 11-W m 10) and 18 mm yr −1 (P m 03-W m 16) were estimated from the peatland to the underlying aquifer, and vertical inflows to the peatland of 116 mm yr −1 were estimated at P m 07-W m 05 (no estimates of vertical flow were available at Cheinu).
Flows in the lower southwestern portion of the peatland were mostly diffuse so outflows from the Misask peatland could not be measured. Outflows were able to be measured at the Cheinu peatland outlet. The average measured flow rate at this outlet was 730 m 3 d −1 (ranging between 371 m 3 d −1 in June 2019 and 1,150 m 3 d −1 in September 2018), which is equivalent to 1812 mm yr −1 (921-2,855 mm yr −1 ) when divided by the peatland area.

Calibrated Models and Simulated Flows
The two models reproduce the measured heads reasonably well, with mean errors close to zero (0.017 and -0.002 m for Misask and Cheinu respectively), low mean absolute errors (0.141 and 0.154 m respectively), and low RMSE values (0.179 and 0.196 m respectively) (Figure 7). The 0.1 m error interval is of  the same order of magnitude as the mean absolute error and includes most of the measured heads.
In the two models, the calibrated recharge values are similar to the range of P net values estimated for the region, and similar for all types of deposits. The calibrated recharge for the Misask peatland is 445 mm yr −1 and is 426 mm yr −1 for the Cheinu peatland ( Table 2). At Misask, calibrated recharge was 414 mm yr −1 for glacial and thin till, 434 mm yr −1 for washed till, and 431 mm yr −1 for fluvioglacial deposits. For the Cheinu model, the calibrated recharge was 450 mm yr −1 for the boulder field area, 414 mm yr −1 for the washed and glacial till, and 416 mm yr −1 for the fluvioglacial deposits.
Calibrated K h values for glacial till and fluvioglacial deposits are within field measured ranges (Table 3). At depths where no measurements were available, calibrated K h values correspond to the expected ranges for the different types of mineral deposits (Domenico and Schwartz, 1998), and show decreasing trends with depth, as expected in unconsolidated sediments. As mentioned above, K h values for the top 1 m of peat were set to measure averages, with a single value for each model layer (not calibrated). Below 1 m, peat K h were set to gradually decrease until 1 × 10 -8 m s −1 at maximum peat thickness.
At the Misask peatland, peatland recharge (447 mm. yr −1 ) dominates inflows (78% of water input), while the groundwater contribution to the peatland (129 mm yr −1 ) represented 22% of all peatland inflows ( Table 2). Simulated groundwater inflows were distributed homogeneously all around the peatland. Results show that seepage of groundwater to the surface of the peatland (490 mm yr −1 ) was the largest outflow from the peatland (85% of water output), while groundwater output through the saturated zone (86 mm yr −1 ) represented 15% of water output. Groundwater outflow was exfiltrated from the peatland relatively homogeneously around its periphery. It is important to highlight here that, to simplify the estimation, vertical and horizontal flows were not distinguished when estimating simulated groundwater inflows and outflows to and from the peatland.
The Cheinu recharge (427 mm yr −1 ) represented 44% of water input to the peatland, while the proportion of inflow from the surrounding aquifer comprised 56% (544 mm yr −1 ). The portion of the peatland that lies on glacial till received 86% of the total groundwater inflow, while the portion of the peatland overlying fluvioglacial deposits received most of the rest of the groundwater inflow. Water outputs were dominated by seepage of groundwater to the surface of the peatland, which accounted for 74% of total outflows (717 mm yr −1 ), while groundwater outflow from the peatland to the aquifer represented 26% of total outflows (254 mm yr −1 ). These outflows were distributed relatively uniformly all around the peatland. Again, horizontal, and vertical flows were not distinguished.
At the Cheinu peatland, the simulated seepage, which is to represent surface flow at the peatland outlet, is markedly lower than the measured flow rates. This discrepancy could be an indication that the outlet stream drains surface and subsurface areas located outside of the peatland watershed in its lower section. Because outflow rates were not measured at the diffuse outlets of the Misask peatland, the simulated outflows could not be validated.

Peat Hydrogeological Properties and Water Table Position
Hydraulic conductivity results from the peat cores (peat K profiles) showed marked differences, despite similar peat characteristics (bulk density, C:N ratios and vegetation assemblages; Robitaille et al., 2021), regional climate, and geological context. At the Misask site, the entire top 1 m of peat showed average values of K h > K v , indicating that lateral flow dominates over vertical flow. At the Cheinu site, the vertical anisotropy resembles that reported by Rosa and Larocque (2008), with K h < K v in deeper layers (>~50 cm), pointing to the importance of vertical flows within the peat, below the acrotelm (40 cm). It is difficult to explain this difference in vertical anisotropy between the two sites considering their similarities, but greater K v in deeper layers at Cheinu could be due to the presence of macropores within the peat enhanced by a larger inflow of mineralized groundwater. This hypothesis is supported by results from Hare et al. (2017), who showed that peatlands with steeper margins (such as Cheinu) are more likely to develop additional preferential discharge areas with macropores and pipe flow, facilitating exchanges through peatland bottom.
As expected, the observed decreasing K with depth and the range of values compared reasonably well with those of other studies (Beckwith et al., 2003;Surridge et al., 2005;Rosa and Larocque, 2008;Bourgault et al., 2018). The decreasing trends with depth are similar to those reported by Beckwith et al. (2003) and Surridge et al. (2005) for peatlands located in the United Kingdom. However, Rosa and Larocque (2008) and Bourgault et al. (2019) found a greater range of K within the first 30 cm of the peat columns in southern Quebec peatlands (i.e., reaching 10 −1 to 10 -4 m s −1 , and sometimes lower K values at 1 m, decreasing by up to six orders of magnitude). The difference between results from this study and those for more southern Quebec peatlands could be due to a dominance of herbaceous peat below the acrotelm at the studied sites (Robitaille et al., 2021).
Relatively stable and shallow WTDs were measured in the peatland sections where peat was the thickest. Waddington (2011) andLukenbach et al. (2015) have already shown that groundwater-dependent peatlands present nearsurface water tables. In this study, the average WTDs of the central wells (18 cm at W m 04 and 15 cm at W c 28) are shallower than the average values of approximately 30 cm reported by Bourgault et al. (2019) for southern Quebec (between 54 cm below ground to 1 cm above ground, with marked temporal variations). The shallower WTDs in the current study could be explained by the short and cool growing seasons that limits evapotranspiration from these peatlands.
For the two peatlands, the results showed more WTD variability at the forested margins where variations in groundwater levels within the aquifer directly influence the peatland water table (Figure 4). This corresponds to the results of Moore et al. (2021) who have shown that areas with shallow peat thickness experience greater water table fluctuations and are susceptible to lose their water table during drought conditions. These marginal areas are frequently characterized by a lagg margin with minerotrophic vegetation (Howie and Meerveld, 2011;Paradis et al., 2015), where groundwater inflow from the aquifer can be substantial (Glaser et al., 1997). Ferlatte et al. (2015) found similar spatialization of connectivity, which peaked at the peatland margin and was attenuated or delayed in the central zone.

Contrasting Water Budgets
The results confirm that the two peatlands are flow-through peatlands, i.e., they receive groundwater from their upstream portions, transmit water through their upper more permeable layers of peat, and discharge groundwater at their downstream edge. However, simulated inflows at Misask were dominated by recharge while inflows at Cheinu had a similar proportion from recharge and from groundwater. The larger proportion of inflows from groundwater at Cheinu is similar to that of Quillet et al. (2017) for slope peatlands. A similar result is also reported by Autio et al. (2020) who simulated larger groundwater flow into peatlands when the slope contrast between the peatland and its immediate uphill surroundings was pronounced. The contrast in inflow sources at the two peatlands could also be due to their position within the watershed, with Cheinu being located close to the regional outlet and thus accumulating more flow from upgradient areas. The larger peatland watershed of Cheinu (0.487 km 2 ) considered to be the immediate groundwater contribution area, as well as the marked topography gradient at the head of the Cheinu peatland (6.7% at Cheinu compared to 0.5% at Misask) likely also contribute to its larger groundwater inflow. The higher proportion of pool area at Cheinu could be a result from its larger catchment area (White and Payette, 2016). Interestingly, simulated seepage and groundwater outflows occupied relatively similar proportions at the two peatlands, although the absolute simulated flows were larger at Cheinu. These results concord with those of Hokanson et al. (2020) who showed that peatlands in the Boreal Plains ecozone of Canada can be potential sources of groundwater to adjacent forestlands. They are also coherent with those of Wells et al. (2017) in the same ecozone where peatlands can generate important runoff.
The Darcy-calculated horizontal groundwater inflows and outflows (eq. (1)) using measured hydraulic conductivities and hydraulic gradients, and either average or maximum peat thickness, captured the simulated groundwater inflows and outflows relatively well. This is an interesting observation since quantifying groundwater inflows and outflows using Darcy calculations can be prone to large uncertainties due to the difficulty in measuring representative hydraulic gradients and hydraulic conductivities (both in the mineral deposits and in the peat). The three estimated vertical flows calculated at Misask (eq. (2)) indicated limited vertical outflow from the peatland to the aquifer in the northwestern and center-northern sections of the peatland, and high vertical inflow in the southeastern sections. Vertical hydraulic gradients of the same order of magnitude and in the same direction as those observed at P m 11-W m 10 and P m 03-W m 16 and similar vertical fluxes are reported by Ferlatte et al. (2015) and Goodbrand et al. (2019). Goodbrand et al. (2019) also report occasional upward hydraulic gradients, but not of a magnitude similar to those observed at P m 07-W m 05. These stations may be located too far apart for a reliable estimate of vertical inflows. Nevertheless, because the peatlands have developed over high permeability mineral materiel, vertical water flows between the peat and the aquifer are probable (Reeve et al., 2000). Vertical flow reversals during sustained dry periods have been reported elsewhere (Fraser et al., 2001;Ferone and Devito, 2004;Ferlatte et al., 2015) but were not systematically observed here.

Connectivity Response to Recharge
The analysis of WTD and flows changes following recharge modifications provides a means to estimate peatland susceptibility to long-term hydrometeorological changes. The simulated WTD changes at the Misask peatland (a decrease of 1 cm to an increase of 3 cm, for recharge changes between −50% and +20%) are very limited, indicating low sensitivity of the peatland hydrology to a changing climate. At Cheinu, the average simulated WTD change is larger (a decrease of 9 cm to an increase of 2 cm), but still relatively limited. These relatively small changes in WTDs reflect the fact that the studied peatlands are in conditions of precipitation excess. Due to the relatively high water table, a drastic drop in recharge would be necessary to lower their water levels sufficiently to accommodate all the direct incoming flow from net precipitation. This study also shows that recharge affects more WTDs in the surrounding aquifer than within the peatlands. This impact is transferred to the peatlands through modifications in hydraulic gradients in the peatland vicinity (upgradient and downgradient) which results in lower (or higher) groundwater inflows and outflows. Autio et al. (2020) have also reported that long-term groundwater drawdown (such as that caused by a recharge reduction) may not be reflected in peatland water levels, due to changes in the proportion of surface and groundwater flows and velocities.
At the Misask peatland, changes in recharge induce a variation in seepage, but a lesser variation in groundwater inflows and outflows. This might result from the relatively smaller watershed and relatively smaller catchment aquifer area of this peatland. In contrast, with its downgradient position and larger catchment area, the Cheinu peatland probably collects more groundwater inflow from deeper portions of the aquifer. These factors probably converge to make the Cheinu peatland inflows and outflows more sensitive to recharge. It has been shown elsewhere (e.g., Winter, 2001) that wetlands dependent almost exclusively on precipitation for their water supply are highly vulnerable to changes in climate conditions. Peatlands that receive substantial volumes of groundwater inflow are less vulnerable because of the large buffering capacity of aquifers to climate change. Hokanson et al. (2018) have confirmed this for drought conditions and shown that flow-through peatlands that intersect the water table are the least vulnerable to deep peat fires. In addition to this, the current work points to the fact that in peatlands with larger groundwater catchment areas, the impact of a change in recharge is also modulated by the impact this change has on the connected aquifer. Peatlands that are strongly connected to aquifers and receive larger fluxes of groundwater as their water fluxes might thus be more vulnerable to changes in recharge conditions, but this still needs to be tested further.
Temperatures and precipitation increases are expected under a changing climate, but precipitation could show more year to year variability and contrasts. Dubois et al. (2022) have identified thresholds of intra-annual precipitation and temperature changes that trigger reductions in recharge for southern Quebec. A similar analysis remains to be done for boreal regions for which it is still difficult to estimate how changes will translate into intra-and inter annual changes. The ecohydrological changes triggered by the wetter and cooler conditions of the Neoglacial period reported by Robitaille et al. (2021) could be quite different in future wetter and warmer conditions.

Peatland Hydrological Functions
It has long been recognized that wetlands are specific components of groundwater flow systems and that their connectivity is affected by their position in the landscape, the aquifer geology, and the climate setting (Winter 1999;Devito et al., 2017;Bourgault et al., 2019). This connectivity can be expressed in terms of hydrological functions of storage, transmission, and runoff (Spence et al., 2011;Goodbrand et al., 2019). Considering WTD and flow sensitivity to recharge and considering the different geomorphic contexts of the two peatlands, it is hypothesized that, on average, the Misask peatland is more dependent on seepage than the Cheinu peatland where the transmission function dominates. Spence et al. (2011) and Goodbrand et al. (2019) have shown that the main hydrological function alternate through the year. Goodbrand et al. (2019) suggested WTD thresholds for these switches between 5 and 15 cm in the Boreal Plains ecozone, i.e., implying that below this threshold peatlands switch from the runoff generation function (corresponding here to seepage) to a groundwater transmission function. In the current study, WTDs are close to these threshold values in the thickest portions of the two peatlands, but larger WTDs were observed in downgradient and marginal portions of the two peatlands. This is an indication that the runoff-to-transmission threshold probably varies within a peatland and between peatlands depending on geomorphic context size of catchment area and local flow directions. The thresholds for peatlands in north-central Quebec might also be different than those of Goodbrand et al. (2019), due to the highly permeable underlying and surrounding geological material which contrasts with the clayey substrate of the Boreal Plains.
It has been argued above that sensitivity to recharge is linked to the peatland geomorphological setting, to the size of the catchment area and to its position in regional flow. These conditions cannot be dissociated from the dominant hydrological function and probably also influence the wetness conditions that trigger switches between functions through the year. For example, at Cheinu, a lowering of the water table could favor peatland transmission at the expense of runoff, reducing the extent of pools at the surface and reducing flow rates at the peatland outlet. In contrast, the runoff function could be more resilient to change at Misask. White and Payette (2016) have shown in fens of northern Quebec that watershed area facilitates pool creation and expansion, and that pools are expected to increase in size and number if climate change brings more humid conditions. Although this needs to be demonstrated further, this work underlines the importance of characterizing the peatland watershed and groundwater catchment area to properly assess peatland hydrological functions, and their sensitivity to hydrometeorological changes.
The dominant type of hydrological functions in peatlands also affects water geochemistry and peatland vegetation (Larocque et al., 2016). Over long periods Belyea (2009) has shown that increasing WTDs can induce major shifts in peatland ecology, from typical ombrotrophic peatland vegetation to a wetherbaceous environment. In contrast, wetter, and cooler conditions such as those evidenced by Robitaille et al. (2021) for the Neoglacial period can trigger ecohydrological changes towards a shift from Sphagnum to sedge-dominated peatlands, and during warmer conditions from herbaceous to Sphagnumdominated peatlands, as evidenced from 33 peat cores from boreal peatlands (Magnan et al., 2022). These long-term variations of hydrological functions can trigger changes in peat accumulation rates (Bridgham et al., 1999). This has been demonstrated at the two studied peatlands where recent warming has increased apparent peat accumulation and induced important paludification (Robitaille et al., 2021). However, feedback mechanisms between hydrology, vegetation, and peat accumulation are not yet entirely understood and their combined modelling remains a challenge.

Model Advantages and Limitations
This work confirms that the UZF package is a useful package to simulate the proportion of overland flow compared to recharge on organic deposits in north-central Quebec. Without this package, it would have been necessary either to use drains all over the peatland to limit water levels or empirical runoff ratios to represent the actual volume of water that reaches the saturated level in the peatlands. This study thus adds to the work of Feinstein et al. (2020) with UZF in Wisconsin fens to demonstrate the usefulness of this package.
Steady-state models have been useful to quantify relatively well average flow conditions within the two peatlands. These models provide a useful macro-scale assessment of the peatland sensitivity to recharge in boreal conditions. It is clear however that transient-state simulations would be necessary to verify hypotheses concerning how the peatland stores, transmits, and releases water throughout the year and how these functions can be expected to change in a changing climate. Because they do not consider transient-state storage, the models used in this study produced linear relations between recharge and WTD and fluxes. However, intra-year variations of storage and discharge through seepage and groundwater flow can vary significantly within peatlands (Oswald et al., 2011;Li et al., 2019), and these changes have reported to be linked to variations of WTD (Kvaerner and Kløve, 2008;Carrer et al., 2015), leading to non-linear relations between recharge, WTD and fluxes. Nijp et al. (2017) have shown that including self-regulation processes through moss water storage and peat volume change can lead to less severe impacts of simulated recharge variations due to climate change. Changing climate conditions are expected to induce modifications in late fall and early spring recharge conditions, and winter recharge events (Jyrkama and Sykes, 2007;Okkonen and Kløve, 2011). Increasing rainfall intensity (Sillmann et al., 2013;Ouranos, 2015) could also modify recharge patterns throughout the year (Fu et al., 2019;Dubois et al., 2022). Jaros et al. (2019) and Autio et al. (2020) argued that aquifer-peatland interactions are best studied using integrated groundwater-surface water flow model. However, this type of model requires an extended data set that was not available in this study, for example storage coefficients for the organic and mineral materials, and continuous transient state surface flow rate data at the outlet.
Finally, some of the simplifications introduced in the models to represent boundary conditions and the heterogeneity of geological and peat materials might have influenced simulated peatland-aquifer exchanges, as well as flow directions and magnitudes between the peat and the aquifer. For example, it has been shown elsewhere that peat hydraulic conductivities at the peatland margins are lower than those in the central portion (Howie and Meerveld, 2011. According to Lapen et al. (2005), this can limit lateral flow to the aquifer and contribute to raising the water table at the center of a bog. Because spatial heterogeneity was not characterized at the two sites (only three peat cores available), it was not included in the models. In the contrasted environment of prairie wetlands, van der Kamp and Hayashi (2009) have shown that deep groundwater can contribute to wetland inflow. Deep groundwater inflows could also be present in the studied peatlands, especially at Cheinu which is located close to an important lake and downgradient from regional flow. However, the absence of measured heads from deeper within the aquifer, and the absence of natural tracer data (e.g., stable isotopes of water) makes this impossible to verify.

CONCLUSION
This research quantified aquifer-peatland hydrological connectivity and identified controlling factors in two pristine boreal peatlands. Considering that the two peatlands have a similar history of peat initiation, are exposed to the same climatic conditions, and have developed over similar glacial deposits, contrasting measured and simulated results emphasized the important role of local factors on their hydrogeological behavior. Measured hydraulic conductivities showed different trends of vertical anisotropy for the two peatlands, but WTD analysis showed stable near-surface levels within the central ombrotrophic areas and deeper water tables in marginal fens area at the two sites.
The groundwater flow models contributed to verify flow hypotheses, quantify unmeasured flows, and explore recharge scenarios leading to changes in connectivity. This work provided indications that, at the Misask peatland, inflows were dominated by recharge and the contribution of groundwater was relatively limited. Outflows to the aquifer were small and mostly exfiltrated by surface seepage. The Cheinu site had a distinct water budget, with inflows dominated by groundwater and most of outflow as seepage. The Misask peatland shows limited sensitivity to changes in recharge, while the Cheinu peatland reacts more markedly to the same changes. It is hypothesized that the greater groundwater dependence of the Cheinu site can be explained by the peatland geomorphological setting, by the size of its catchment area and by its position in regional flow. These factors necessarily influence the hydrological functions of the two peatlands and their sensitivity to long-term recharge changes. This work underlined the importance of the peatland watershed and groundwater catchment area to assess its resilience to hydrometeorological changes.
More work needs to be done to confirm results from this study, notably the influence of transient processes which could trigger feedback processes leading to either stable long-term conditions or in contrast drastic changes in hydrological conditions leading to changes in vegetation and peat accumulation. To achieve this, long-term transient-state simulations including flow, vegetation and peat accumulation need to be performed. This work has shown that aquifer-peatland connectivity at the regional scale needs to be included in this assessment.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
CL, ML, and SG contributed to develop and implement the field and modelling approach. ML supervised the project and MG contributed for the ecohydrology components of the project. Data management, simulations, and figure preparation were done by CL. MG and ML obtained the research grant, All authors contributed to write the paper.