Quantifying groundwater fluxes from an aapa mire to a riverside esker formation

Diverse peatland–aquifer interactions can be expected in boreal mosaic quaternary landscapes. In recent years, there have been major efforts to identify temporal and spatial connections between aquifers and peatland (Ferlatte et al. 2015; Quillet et al. 2017). Groundwater (GW) is generally an important part of the water balance of fens (Boeye & Verheyen 1992; Feinstein et al. 2019). Local geomorphological and hydrogeological settings (Lowry et al. 2007; Kværner & Kløve 2008; Isokangas et al. 2019) influence runoff generation and spatial variations in seepage within peatland complexes (Hare et al. 2017; Isokangas et al. 2017). At peatland margins or at the boundary with mineral soil formations, peatlands can also supply water to underlying mineral deposits (Reeve et al. 2000; Ferlatte et al. 2015). In recent years, several detailed field measurement campaigns (Lowry et al. 2007, 2009; Isokangas et al. 2019) and detailed simulations of regional GW patterns (Åberg et al. 2019; Feinstein et al. 2019) have increased understanding of these complex systems. However, only a few studies (e.g., Bourgault et al. 2014) have monitored and modeled water flow from peatlands and its connections to regional flow systems.

Peatlands play an important role in the hydrological and hydrogeological dynamics of boreal landscapes (Waddington et al. 2015), influencing flow systems and modifying the ecology and biodiversity of these landscapes (Chapman et al. 2003). Aquifer–peatland connections are typically influenced by peatland location (Lowry et al. 2009) and peat development stage. However, site-specific conditions may have a strong influence, making it difficult to extrapolate general findings to a specific area. There is a strong need to identify factors with site-specific influence and to develop tools for investigating aquifer–peatland connections, particularly in light of increasing land-use pressure in boreal landscapes, in which peatland–aquifer complexes are, e.g., typical drinking water sources (Johansen et al. 2011; Rantala et al. 2017; Xu et al. 2018) or are used for forestry (Rossi et al. 2014) and agriculture (Paavilainen & Päivänen 1995). Unknown pathways and changing connections, especially at peatland margins, might have momentous and undesirable influences on local surface water (SW) and GW flow paths.

This study was conducted to investigate the connection between an aapa mire margin and riverside esker formation. The main hypothesis tested was that in this complex setting, recharge patterns are the reverse of the usual flow pattern from peatland to aquifer reported in the scientific literature for typical peatlands for which hydrology has been studied. The geological and hydrological conditions that drive peatland–aquifer connections are unclear and, in particular, riverside aquifer formations in plain landscapes have received little attention in previous research. These geological formations may play a significant role in local biogeochemical processes and greenhouse gas (GHG) fluxes. Any organic matter (OM)-rich peat water entering a sand formation would alter the GW chemistry, promote biogeochemical processing of OM in GW layers, and eventually affect the local GHG fluxes at the surface. The overall aim of this study was to determine how hydrogeological conditions influence peatland–aquifer connections in an aapa mire complex. Specific objectives were to identify: (i) how the aapa mire margin is hydrogeologically connected to the riverside aquifer and (ii) spatial and temporal variations in recharge of peatland water to GW.

The Tähtelä site (Figure 1) is located in northern Finland, close to the municipality of Sodankylä. Long-term (1981–2010) mean annual temperature in the area is −0.4 °C, and mean annual precipitation is 527 mm (Finnish Meteorological Institute, FMI). Nearly half of annual precipitation falls as snow. Snowmelt typically occurs during the second half of May, and the first snow appears in mid-October. Long-term mean annual evaporation (1961–2010) of 315 mm has been measured by an evaporation pan (Moroizumi et al. 2014). The landscape is rather flat and is dominated by aapa mire complexes, while the geology next to the river Kitinen is dominated by sandy eskers and underlying bedrock or glacial till deposits. The altitude of the mire is about 180 m above mean sea level (MAMSL), and hills in the area rise to 5–10 m above the surface of the mire. The area is located close to the Central Lapland ice-divide zone, which means that surficial deposits are thicker and more complex compared with the average in Finland (Johansson & Kujansuu 2005). The peatland margin areas consist of a paludified zone of mineral soils on the edges of the aapa mire, and the aapa mire next to study areas is treeless but contains few poorly growing trees. The mires can be considered ‘pristine’, although some old ditches exist in the area. South-east part of the peatland is a poorly drained peatland forestry area.

In total, 11 GW wells (GWWs) were installed at Tähtelä in 2013 to observe the GW flow process. Two additional GWWs were installed at the nearby peatland in 2014. In 2015, the water level during summer (June–September) was monitored manually in all wells at 2-week intervals, and continuously (15-min intervals) with water-level loggers (Solinst levelogger) installed in four wells (GWW0, GWW2, GWW4, and GWW9). GWWs were made from galvanized steel (inner diameter: 8 cm) and screened for the first 1 m from the tip. GWWs were designed to follow the hypothetical GW flow path (Figure 1) from aapa mire to river and, thus, measure changes in local GW levels. GW flow path assumption was based on pre-analysis of the digital elevation model (DEM, 2 × 2 m) and field measurements. GWW0 was 2.5 m deep, GWWs 1–10 were 4.5 m deep, GWW11 was installed to mineral layer under the peat to 2.5 m depth (peat depth: 1.3 m), and GWW12 was installed only on peat layers at 1.5 m. GWW6 was not used to monitor GW levels or for sampling. Meteorological conditions were measured at a nearby official FMI weather station.

The hydrogeological structure of the area was studied using a variety of methods. Soil type was investigated during GWW installation, and it was found that the esker contained a layer of homogeneous fine sand that was several meters thick. These measurements were supplemented with additional manual peat drillings and ground-penetrating radar (GPR) measurements (Malå 100 MHz system) in 2016 and 2017. The total length for GPR lines was approximately 15 km, and the MALÅ Object Mapper software was used to post-processing the GPR lines. The spatial distribution of peat thickness and underlying geology were interpolated for the area using the natural neighbor method with GPR measurements.

All GWWs (except GWW6, based on previous information) were sampled for stable isotopes of water (δ²H, δ¹⁸O), dissolved organic carbon (DOC), and dissolved inorganic carbon (DIC) at 2-week intervals during summer 2014–2017 from May to October. All samples were pumped from GWWs, by filling 50-mL sterile plastic bottles with no headspace, and first pumped water was not taken to the sampling bottle to avoid ‘standing water’ in the pipes. In addition, snow core samples were taken during spring winter from the peatland area (next to GWW10) to determine isotope values. Precipitation samples were also collected at the official Global Network of Isotopes in Precipitation (GNIP) sampling points in Rovaniemi and Sodankylä. All samples were stored in darkness at +4 °C until analysis. Dual isotope ratios (²H/¹H and ¹⁸O/¹⁶O) in water samples were determined using cavity ring-down spectroscopy with a Picarro L1102-i spectrometer at the University of Oulu, Finland, and calibrated using in-house standards calibrated to Vienna Standard Mean Ocean Water (VSMOW). All isotope ratios are expressed in δ notation relative to VSMOW, with precision for δ¹⁸O and δ²H values of ±0.1 and ±1.0‰, respectively. The DOC and DIC samples following standard SFS-EN 1484:1997 were analyzed at the accredited laboratory of Natural Resources Institute Finland (LUKE).

The GW flow from the aapa mire complex to the esker was simulated with the MODFLOW GW model (McDonald & Harbaugh 1988). It consists of a one-layer model (no continuous stratigraphy was detected in the esker) with a uniform cell size of 20 × 20 m² and a total model area of 0.53 km². The main purpose of the model was to estimate the total water flow volumes from peatland to the sandy esker. The single model layer was set to the model, since GPR analysis revealed that esker soils did not have any layer between the soil surface and bottom material. Boundary conditions fell into two categories. Constant head boundary conditions were applied to the eastern model boundary along the peat and to the western model boundary along the river. Based on results from GPR measurements, the eastern model boundary was set along the location where hydrological exchange between peat and esker was assumed and measured. The study assumed a constant head for the peatland and margin area. This assumption was supported by measured water table data from peatland GWWs, which showed minimal variation in the water table. In the model, we set the model boundary to GWW10 which represents the peatland margin area and thus at the same time the water table in peatland, but also the water table input to the esker system. Other boundaries were considered no-flow boundaries, and it was assumed no major flow paths from the sides. This assumption is supported by observed water levels and gradients inside GWW5 and GWW7, which were similar to those in other GWWs, suggesting no side boundary GW flow. Along the river, the constant head for the whole thickness (simple 1-layer model) was assigned a value of 174.6 m, derived from the elevation at the river surface in the DEM. The water level in river Kitinen at the site location is fairly constant due to damming of the downstream hydropower plant. Surface elevations were derived from a DEM (2 × 2 m²), and bottom elevation followed permeable layer topography measured with GPR. For the geological structure, one homogeneous layer of sand was assumed to cover the whole modeled area, with impermeable bedrock (or till) underneath. The bottom elevations of this single layer (and therefore the bottom elevation of the 3D model) were obtained by the interpolation of data from GPR measurements of the area (see ‘Field measurements and data collection’). The parameter estimation program PEST (Doherty 2019) and manual adjustment were used to adjust the hydraulic conductivity parameter.

Steady-state and transient models were tested. In the steady-state model, the constant head along the peat was assigned the average GW level measured in GWW10 (177.56 m). In the transient model, the GW levels recorded in GWW10 were assigned to each period as a constant head along the peat. Continuous GW measurements from wells GWW0, GWW2, GWW4, and GWW9 were used to calibrate the models. For the steady-state model, the average piezometer levels were used as the calibration target. Recharge for the steady-state model was assumed to be 0.00075 m d⁻¹ (based on the average annual precipitation of 520.5 mm at the FMI station), and actual annual evapotranspiration was assumed to be 248 mm for the area (estimated from regional weather stations).

The steady-state model was converted into a transient model by applying the following changes: 10 periods matching the manual GW measurement intervals, covering a total of 453 days, were defined. These stress periods were subdivided so that the model was run in time steps of approximately 1 week. The first period (1 September 2013–11 March 2014) followed the settings for the steady-state model with monthly time step and was only used to adjust the model before actually applying time-variant observed data. For the transient model (2 weeks' time step), local daily meteorological data from the FMI station were used to estimate recharge variations over time using a simplified water balance approach. Values for snow water equivalent (SWE) were taken from the two closest points (Sodankylä Unari and Raudanjoki, about 30 km away) for which measured data were available from Finnish Environment Institute (SYKE). Snow accumulation and melting were added as changing total recharge for each time step based on measured SWE, where infiltration rate was estimated from ET.

The GPR measurements revealed a geological subsurface structure that allowed the recharge of DOC-rich peatland soil water to underlying sand layers and to the riverside esker formation. Side areas of the peatland consisted of organic peat layers up to 1 m thick above the sand formation. In the assumed GW flow direction, there was a clear depression in the bedrock or glacial till that formed a sand pipe toward the river (Figure 2). Sand layer depth increased toward the river, creating the riverside esker and conducting GW to the river Kitinen. Below the sand layer, bedrock or hard till layers formed an impermeable boundary. Topography and geological formations are known to control SW and GW patterns, and thus the formation of the water-logged environments required for peatland development and survival. In aapa mire complexes, peat type can vary spatially from ombrotrophic to minerotrophic, but flow patterns in fen types are typically from higher GW head in mineral soil toward lower peatland formations. The present study revealed another type of connection, where an aapa mire complex margin is connected to a riverside esker formation through a local subsurface sand formation.

The GW flow direction within the study area was toward the Kitinen river and the GW level varied between 177 and 175 MAMSL (Figure 1). The GW fluctuation in the esker formation was rather stable and was confirmed by MODFLOW simulations. In the GW model, the calibration targets for four wells (GWW0, GWW2, GWW4, and GWW9), which were used in the model, were replicated accurately in the steady-state model (mean residual: 0.12; mean absolute residual: 0.26; root-mean-squared residual: 0.35; mean squared error (MSE): 0.12). The parameter estimation program PEST (Doherty 2019) and manual adjustment were used to adjust the hydraulic conductivity parameter for a Kf value of 8.9 × 10^–4 m s⁻¹, which is typical for sandy soils. The calibrated steady-state model was then run as a transient model (mean residual: 0.22; mean absolute residual: 0.38; root-mean-squared residual: 0.56; MSE: 0.31). Temporal variations in GW levels were accurately replicated with measured GW values by the transient model (Figure 3 and 5(b)). As modeled results showed good agreement with measurements, no further calibration was done with the transient model. In the transient model when compared with measured values, GWW0 showed the highest uncertainty (Figure 3). The transient model showed the highest deviations during the snowmelt period, at least for GWW0 where water-level fluctuation at nearby bond is also influencing on the GW levels. We used a simple model approach with constant head boundary conditions, single layers both in the peat and in the sand, and a simplified approach to estimate recharge. Using a more complex model could have resulted in more accurate modeling than our simplified approach. However, it is important to note that our simplified model produced good estimates from seasonal fluctuations and it was thus further possible to estimate GW seepage fluxes from mire to esker area.

The flow model confirmed that the major flow direction from the peatland is toward the river Kitinen (Figure 4), as also seen in the GW head observations. According to the model, the main GW recharge period was during snowmelt, and there was only minor recharge from summer precipitation (Figure 5). In the study area (0.53 km²), water flow from peatland to sand esker was estimated to be 873 m³ d⁻¹, or on average 70% of daily GW outflow in the esker area, thus representing a major part of esker water budget. The transit time using particle tracking in the GW movement from peatland recharge to the river Kitinen based on modeling results was estimated to be 692 days. Annual variation in the GW level in the riverside esker followed a typical pattern for boreal conditions (Mäkinen 2003; Åberg et al. 2019), with a maximum in June after peak snowmelt and a minimum at the end of winter (April). Recharge values in the Central Lapland have been reported to vary from 42–61% of annual precipitation in a sandy esker (Hyyppä 1962) to 3–41% of annual precipitation in a mire located around 20 km north of the Tähtelä site and representing rather similar conditions (Åberg et al. 2019). At this study site, recharge represented on average 30% of annual precipitation. GW fluxes vary notably depending on local geology and size of GW formation, and thus fluxes from mire to esker cannot be directly compared with other sites. However, GW heads showed similar patterns that have been observed in previous studies in northern Finland (Okkonen & Kløve 2011; Rossi et al. 2012).

The Tähtelä site is an excellent example contradicting the oversimplified assumption that the recharge rate in eskers is only directly proportional to precipitation. Future studies seeking to determine recharge rates in similar aapa mire complexes should thus include the analysis of the local geological structure, especially for the case of interfaces between upstream aapa mires and mineral soils.

Measured δ¹⁸O values in the mire well (GWW10) were close to those in summer precipitation (mean summer precipitation −11‰) but became more GW-dominated from GWW8 onwards the river Kitinen (Figure 6). Both GWW10 (in peatland) and GWW9 (in mineral soil next to peatland) showed isotopic signals of peat water (summer precipitation) and evaporation (d-excess) (Figure 6). These signals disappeared on approaching GWW5 and GWW4 in the middle of the esker system. River isotope signals differed notably from the closest GWW0 (Figure 6), indicating no recharge from river to riverside esker. Mean measured δ¹⁸O and δ²H values in GWW0 and GWW8 (−14.18 and −104.98 ‰, respectively) were close to values reported previously for regional GW (mean: −14.68 and −107.8‰, respectively) in Siurunmaa, Sodankylä (Kortelainen 2007), whereas those in GWW9 and GWW10 (mean: −11.86 and −90.48‰, respectively) differed slightly from values reported for open mire sites in Viiankiaapa, Sodankylä (−12.28 and −94.4‰, respectively) (Åberg et al. 2019). There were clear seasonal variations in δ¹⁸O values in wells closer to peatlands, where snowmelt inputs gave strongly negative values (−17‰ δ¹⁸O), but these shifted to follow summer precipitation values already in June. The seasonal variation was not as strong in wells closer to the river Kitinen, but inputs from snowmelt still had a discernible effect (Figure 6). In these, the GWW water level was more stabilized from the GW system. Similarly, small changes in summer isotope values were apparent in lower GWW, but continuous GW level monitoring indicated no major summer precipitation recharge directly from sand to esker. The GW level declined gradually during summer, and thus the summer δ¹⁸O signals in GW most likely originated from recharging peat water to the GW system (Figure 7), rather than from precipitation or regional GW patterns. At peat margin and close by areas, GWW (GWW10 and GWW9) isotope values and d-excess showed (Figure 7) evaporation signal (values below local meteorological water line, LMWL) and supported the hypothesis from recharging of water from peatland to esker through margin area.

Eskers are natural filter systems affecting the physical and biogeochemical processes in GWs (Lindroos et al. 2002). Analysis showed the transformation of DOC in GW to DIC along the rather short travel path from GWW10 to GWW7, most probably due to bacterial metabolic processes (Stegen et al. 2012). The DOC and DIC concentrations varied markedly along the GW gradient from the peatland to the river Kitinen (Figure 6(c) and 6(d)). Both GWW10 (up to 70 mg L⁻¹) and GWW9 (mean: 4.94 mg L⁻¹) showed high DOC values, while DOC values were lower at most other GWWs (mean: 1.17 mg L⁻¹). GWW11 was in sand below the peatlayer, and GWW12 was only in peat. DOC values in GWW11 were at the same level as in GWW9 indicating minor water transfer between peat layers and underlying sand. However, at peat margin GWW10 DOC values increased notably indicating water movement and DOC transfer from peat layers to sand deposits. The DIC values increased from GWW10 to GWW9, and then gradually decreased toward the river, indicating the transformation of DOC to DIC and of DIC to carbon dioxide (CO₂). Another explanation would be the dilution of DOC values and/or the transport of DIC from side boundaries to the well, but this hypothesis is not supported by water-level observation or modeling. Water levels and gradients inside GWW5 and GWW7 were similar to those in most of the other GWWs, suggesting no side boundary GW flow. This indication is supported by the observed positive annual net ecosystem exchange of CO₂ above the forest, suggesting that the forest site is a net source of CO₂ despite the growing tree stand (e.g., Owen et al. 2007). The DIC concentrations were highest at GWW9, GWW8, and GWW4 and decreased already in GWW2. However, wells located to the sides of the main GW transect (GWW5 and GWW7) had lower DIC concentrations, indicating that carbon transport from the peatland followed a restricted GW pathway (Figure 1).

The present study is one of few studies to demonstrate the flow pattern from mire to esker (Hokanson et al. 2020). Typically, eskers are assumed to be water sources for peatland mires, but here we show opposing flow patterns. The objective of this research was to better understand how hydrogeological conditions in the aapa mire margin influence peatland–esker recharge patterns. Detailed field measurements, including continuous water table level, water stable isotopes, dissolved carbon, and GW flow model, were used to investigate local water movements. Our main hypothesis on the recharge of DOC-rich waters from the aapa mire complex to riverside esker was confirmed by the data obtained. First, GPR measurements revealed a local sand pipe at the study location, which physically enables flow from peatland to sandy esker. Second, stable isotope and DOC analyses showed peaty water transfer from peatland to esker, with DOC transformed to DIC along the pathway to the river Kitinen. Third, the MODFLOW model and GW measurements were used to quantify flow direction and volume between peatland and sandy esker.

A conceptual hydrogeological model of the study area was formed based on these results (Figure 8). Results and a conceptual model show the importance of considering peat margin areas and local geological structures in GW analysis in aapa mire complexes. Similar SW and GW connection need to be taken into account in the water and land management of boreal eskers.

This study is limited to a small area along the long river corridor, but similar aapa mire conditions are found throughout the boreal and sub-arctic regions. This study provides evidence on the importance to consider peat margin areas and local geological structures in GW analysis in aapa mire complexes and contributed to increasing the understanding of how aapa mire peatlands and eskers are connected. Results of this study point to the importance of considering the hydrogeological properties in water and land management planning and implementation.

This study was funded and supported by the Kvantum Institute Spearhead Project and the K.H. Renlund Foundation.

Data cannot be made publicly available; readers should contact the corresponding author for details.