The reference peatland for simulating groundwater flow and discharge in this study is based on a coastal peatland near the city of Rostock, northeast Germany, the Hütelmoor Nature Reserve (“Naturschutzgebiet Heiligensee und Hütelmoor”; 54.2139 N 12.1725 E) (Figure 1). This coastal peatland borders the Baltic Sea along a 3-km long shoreline, but is physically separated by a 40-m wide dune dike that thins out in its northern part. The area was shaped by the Weichselian ice age. The glacial till is covered by glacio-fluvial sands with 2.5–15 m thickness, which now form a shallow aquifer (Ibenthal, 2020 and references therein). During the Littorina transgression 8,000–1,200 BP, the rising groundwater level resulted in the development of a paludification fen from 7,000 BP onwards (Kreuzburg et al., 2018). Peat thickness is up to 3 m in the central peatland and near the coast, while it thins out towards the forest. A special feature is the outcropping of the peat at the coastline and in the shallow Baltic Sea. Strongly decomposed sedges and reed make up the peat. Elevation of the peatland is −0.15–0.75 masl (meters above sea level) with a total area of 350 ha (Ibenthal, 2020). Annual precipitation, annual evapotranspiration, and average daily temperature are 693 mm, 604 mm (1951–2010; Miegel et al., 2016), and 9.7°C (at Warnemünde, ∼7 km away, 1990–2019; DWD(Deustcher Wetterdienst), 2020), respectively.

The Hütelmoor, similar to many coastal peatlands in the southern Baltic Sea region, has been artificially drained, diked, and utilized as pasture. Initial drainage started in the 18th century and extensive pumping from 1976 to 1991 led to degradation and compaction of the peat’s upper decimeters (Jurasinski et al., 2018). In December 2009, the construction of a ground sill with an elevation of 0.4 masl at the outlet of the peatland’s ditch system enabled rewetting and renaturation of the peatland. Since then, most of the peatland is flooded with surface waters. The ditches still drain the peatland when the water level exceeds the groundsill, they affect local groundwater flow fields and divert water inflow from the forest (Ibenthal, 2020). The peatland is still separated from the Baltic Sea by a dune dike, that is not maintained anymore. The electrical conductivity (EC) in ground- and surface water in the peatland ranged from 4–13 mS/cm and thus revealed brackish waters, originating from earlier inundations and occasional seawater inflow via the dike system during storms (Ibenthal, 2020). In comparison, several EC measurements from the nearshore waters in front of the Hütelmoor ranged from 18 to 27 mS/cm, corresponding to a salinity of 11–15 psu. Due to the brackish nature of groundwater in the peatland, we name the submarine groundwater discharging from land “terrestrial” (and not fresh) SGD in this study.

We simulated terrestrial submarine groundwater discharge in the Hütelmoor using the HYDRUS-2D modeling package. HYDRUS simulates water movement in variably saturated media by solving the Richards equation for Darcian water flow: ∂ θ ∂ t = ∂ ∂ x i [ K ( K i j A ∂ h ∂ x j + K i z A ) ] − s (1) where θ is the volumetric water content, h is pressure head, x _i are spatial coordinates, t is time, S is the sink term, K i z A are components of anisotropy tensor K i j A , and K is unsaturated hydraulic conductivity (Šimůnek et al., 2018). The equations are solved using the Galerkin-type linear finite elements method which can be discretized spatially and temporally. An iterative process is used to solve the equations due to the non-linear nature of the Richards equation (Šimůnek et al., 2018).

HYDRUS can simulate the flow of water both in the unsaturated dune dike and the fully saturated peatland. The reference scenario comprises a 2D cross-section reaching from the central part of the Hütelmoor Nature Reserve peatland into the shallow Baltic sea (Figures 1, 2). The topography, geological stratification, material properties, and boundary conditions were based on measured data, and steady-state groundwater fluxes were simulated to yield a long-term average of terrestrial SGD. Fluxes originating from aquifer, peat, and dune were quantified separately. In subsequent model runs, these factors were changed one by one based on actual and realistic predicted conditions, to determine their respective impacts on SGD.

Water inflows (Figure 2–blue arrows) into the coastal peatland system come from 1) the ponded surface waters above the peat, 2) precipitation at dune dike and beach, 3) recharge from recirculated SGD at the seafloor (corresponding to Constant Pressure Head 2), and 4) lateral groundwater inflow from the landside. Outflows (orange arrows) from the system occur as evapotranspiration in the dune and the total SGD from the seafloor. Uncertainty in water flow from the sea-submerged aquifer sands (left boundary) makes it either an inflow or outflow region. Infiltration from ponded peatland surface waters is also expected into the dune dike.

We established the modeling domain (Supplementary Figure 1) using published and available data. It extends from the shallow Baltic Sea to the central peatland (Figure 1). Elevation measurements at the groundwater observation wells (Ibenthal, 2020) and seawater level were used as reference points of the domain, while ground heights in between these reference points were obtained from a Digital Elevation Map (provided by Landesamt für innere Verwaltung Mecklenburg-Vorpommern, Schwerin, Germany) using the 3D Analyst Toolbox of ARCMap. The materials were distributed following the geological profile constructed based on sediment cores (Jurasinski et al., 2018; Ibenthal, 2020). Bathymetric data from marine surveys (Kreuzburg et al., 2018) completed the model geometry, taking note of the longshore bar mound 100 m from the coast. The modeling domain’s total length is 930 m, of which 245 m extend into the sea. In the vertical direction, it extends from 3.7 masl at the top of the dune dike to −9 masl at the aquifer bottom. The bottom boundary follows the surface of the glacial till.

The modeling domain’s discretization was carefully scrutinized and yielded a final target size of 1.0 m, x-direction stretching factor of 25, and a smoothing factor of 1.8. To account for the highly non-linear water flow in the unsaturated zone of the dune dike, a surface mesh refinement of 0.05 m was applied in the upper beach/dune sands. Surface refinements of 0.2 and 0.5 m were also applied at the lower beach/dune sands and peat materials, respectively. A total of 3,454 nodes resulted from this configuration with 351 nodes and 351 elements on the domain boundary.

The maximum model runtime was set at 360 days. All simulations achieved steady-state in the given runtime except for the minimum beach/dune K_s scenario, where runtime was increased to 720 days. Initial, minimum, and maximum time steps were set at 0.01, 0.001, and 5 days, respectively.

We applied the single-porosity van Genuchten-Mualem model. Soil hydraulic properties of the reference scenario (Table 1) were derived from soil data gathered previously at the study site. For soil texture of the aquifer sands and beach/dune sands, average particle size distributions were calculated from 34 to 6 soil samples, respectively. The resulting composition of aquifer sands was 89.4% sand, 7.5% silt and 3.1% clay. The beach/dune sands were sandier with 98.3% sand, 0.9% silt and 0.8% clay. Soil hydraulic parameters were then generated using the neural network prediction “Rosetta” incorporated into the Hydrus software (Schaap et al., 2001). Peat hydraulic parameters were calculated with the pedotransfer functions for sedge bulk density of ≤0.2 g cm⁻³ proposed by Liu and Lennartz (2019), based on measured bulk density ranging from 0.12 to 0.15 g cm⁻³ in the peatland (L. Gosch, personal communication, August 23, 2019) and 0.17 g cm⁻³ in the peat exposed at the coastline (Gosch et al., 2019). The residual water content for peat was assigned a value of 0.2. For the tortuosity and pore-connectivity parameter (L), a value of 0.5 is acceptable for degraded fen peat under wet or saturation conditions (Dettmann et al., 2014; Liu and Lennartz, 2019) as well as for sands, as this is the average value for most soils based on calculations by Mualem (1976) (Simunek et al., 2018). For the peat and aquifer sands, saturated hydraulic conductivities (K_s), values determined with slug tests at the study site were used (Ibenthal, 2020); the resulting K_s values (geometric means) are 8.64 × 10^–3 m d⁻¹ for the peat (n = 4) and 1.73 m d⁻¹ for the aquifer sands (n = 4). Dune K_s of 11.8 m d⁻¹ (n = 8) was laboratory-determined on samples taken in the northern part of the dune and was adapted from a previous study (Mohawesh et al., 2017).

The model boundary conditions (BC) were set as follows (see also Figure 2): Atmospheric BC was applied to the dune surface wherein long-term averages of meteorological variables were used (average precipitation: 1.90 mm d⁻¹; average evapotranspiration: 1.65 mm d⁻¹) (Rostock-Warnemünde 1951–2010; Miegel et al., 2016). At the beach, a seepage face BC was applied which switches to atmospheric BC in the absence of seepage. However, no seepage was observed across the seepage face in any of the simulations. In HYDRUS, “seepage face” refers to the boundary condition applied in surfaces where water leaves the saturated part of the flow domain. Constant head BCs in hydrostatic equilibrium were set on the left side (Constant Pressure Head 1) and at the seafloor (Constant Pressure Head 2), the pressure head corresponds to the average sea water level at Warnemünde, Rostock (2005–2015) of 0.091 masl (WSV(Wasserstraßen und Schiffahrtsverwaltung des Bundes), 2020). Likewise, the right-side boundary (Constant Pressure Head 3), which continues to the peatland surface, was assigned a constant BC in hydrostatic equilibrium, corresponding to the average peatland water level of 0.357 masl (September 2016 to October 2018; derived from 73,735 individual 15-min interval measurements; converted to equivalent freshwater heads). No flux was assumed at the bottom as the glacial till was considered impermeable.

Initial pressure heads were set to correspond with the constant head boundary conditions. As such, the initial bottom pressure head amounted to 9.357 m. A slight slope (0.05, x-direction) from the sea up to behind the dike was added because of pressure head differences. The peatland’s ponded conditions have a higher pressure head compared to the seawater on the left side. The angular slope was not applied in sea level scenarios as the hydraulic gradient’s difference was only 0.001 m.

Most groundwater bodies are heterogeneous with hydraulic conductivity variability influencing groundwater transport (Peña-Haro et al., 2011). In peat soils, it has been shown that a heterogeneous hydraulic conductivity distribution can lead to complex groundwater flow patterns (Beckwith et al., 2003). Since our simulations are based on homogenous hydraulic conductivities in the three materials considered, uncertainty of the model outcome resulting from the variability of K_s cannot be assessed. To address this, we performed 100 extra simulation runs with random K_s values (Supplementary Table 1). The number of replicates of the K_s measurements was too small to derive distributions (aquifer, peat), or samples were only taken from a single site (beach/dune K_s). We therefore additionally analyzed K_s values calculated from particle size distributions following Beyer (1964) both for the aquifer (n = 34) and the beach/dune sands (n = 6). The geometric means were close to the ones of the measured K_s values. All K_s values were then log₁₀-transformed. Random numbers were drawn from normal distributions with means based on the measured values, and the standard deviations of the Beyer-derived values.

For the peat K_s, the coefficient of variation reported by Baird (1997) for a degraded peat was used to calculate the SD. Other published datasets for fen peat (Liu et al., 2016; Wang et al., 2021) were not considered because of much higher K_s values. The log₁₀ mean and SD were then used to generate the random K_s values. In Supplementary Table 1, a descriptive statistics summary of the random K_s values is listed.

To investigate the factors that affect the magnitude and location of terrestrial SGD from coastal peatlands, three clusters of scenarios were analyzed: 1) soil physical properties, 2) hydraulic gradients and, 3) geological stratification and topography (Table 2). The parameters were changed one at a time to better determine their individual impacts on SGD.

HYDRUS calculates flow velocities for each of the nodes and reports integrated fluxes across boundaries. Terrestrial SGD was assumed to be equal to the flux across the Constant Pressure Head two boundary, which extends across the entire seafloor (Figure 2). The net flux of recirculated seawater across this boundary is zero under steady-state conditions. The terrestrial SGD comprises all the water from the peat, aquifer sands, and beach/dune sands that come out as submarine groundwater discharge. Water fluxes from each of the different materials were determined using mesh lines in the material boundaries (Figure 2); fluxes across mesh lines are calculated by HYDRUS in the same way as those across the boundary conditions.

Discharge and recharge regions in the seafloor were identified by looking at velocity vectors and streamlines, both of which are outcomes of HYDRUS graphic user interface post-processes. After identification of discharge and recharge regions, we added meshlines to these regions to quantify water fluxes separately. Discharge determined from these regions can be defined as “total SGD”, since it includes both terrestrial and recirculated SGD. However, caution should be taken when interpreting total and recirculated SGD as the model does not consider marine forcings such as waves. The calculated fluxes across boundaries or mesh lines are reported in m² d⁻¹. To determine the seepage rate from the seafloor (m d⁻¹), the fluxes were divided by the total length of the respective boundary or mesh line.

The upscaled total terrestrial SGD for the entire 3-km coastline of the Hütelmoor is 240 m³ d⁻¹. Lower and upper limits of estimates based on K_s uncertainty are 114 and 596 m³ d⁻¹, respectively. The results are slightly higher than previous terrestrial SGD estimates from the same study site with 15–164 m³ d⁻¹ (Miegel et al., 2016) and 180 m³ d⁻¹ (Ibenthal, 2020), which were calculated using water balance equations of the contribution of peat and aquifer sands (Miegel et al., 2016) and MODFLOW 3D numerical simulations including the dune sands (Ibenthal, 2020).

Comparing with other wetland environments, our average seepage rate of 0.33 cm d⁻¹ is comparable in magnitude to a mudflat in China (Qu et al., 2017) and lower estimates of a tide-dominated coastal wetland in Taiwan (Hsu et al., 2020) (Table 4). However, it is expectedly lower than seepage rates from crab burrows-influenced tropical mangrove systems in Australia due to preferential flow pathways (Tait et al., 2017) and from a subtropical estuarine creek adjacent to a dune-wetland system (Sadat-Noori et al., 2015). The seepage rate near the shore (1.05 cm d⁻¹) is comparable to rates observed in sandy beaches (Kotwicki et al., 2014; Qu et al., 2017) and a sandy pockmark (Virtasalo et al., 2019) in the Baltic Sea. This result may be explained by the fact that terrestrial SGD could also be generated from small dune dike systems with additional inputs from infiltrating peatland surface waters. However, it should be noted that our average seepage rate represents terrestrial SGD only while seepage rate near the shore has contributions from recirculated SGD. Fresh SGD contributes 0.01–10% of surface water runoff (Church, 1996), with recent global estimates of 0.6% (Luijendijk et al., 2020). In tide-dominated systems, a conservative estimate of fresh SGD contribution is ∼5% of the total flux (Li et al., 2009; Santos et al., 2009:; Hsu et al., 2020). Total SGD flux estimates could thus be larger with possible large inputs from recirculated SGD.

This is the first time that HYDRUS-2D was used to determine SGD from a coastal peatland system, and it was able to address our objectives. The model has performed reasonably well—all scenarios have less than 1% relative mass balance error, the prescribed rate for water flow simulations. The lowest and highest relative errors were 0.01 and 0.89%, respectively.

A drawback of the model for applications at the coast is that it does not account for density-driven flow. We assume, however, that density effects can be neglected at our study site due to small density difference between the Baltic Sea and the peatland’s groundwater. While the salinity of the Baltic Sea is comparably low near the study site with 11.4 psu (IOW(Leibniz Institute for Baltic Sea Research Warnemünde), 2020), high electrical conductivities and chloride concentrations have been recorded in the groundwater of the peatland (Ibenthal, 2020). They are rather homogeneously distributed in the peat and the underlying aquifer and are attributed to former floodings of the peatland with seawater, as has been proven with sulfur isotopes (Koebsch et al., 2019). The groundwater in the aquifer of the considered transect is decades old, as has been revealed by tritium-helium dating (Ibenthal, 2020). The densities of water samples collected between 2019 and 2020 from the groundwater peatland (n = 32) and dune (n = 3), from surface water (n = 15) and Baltic Sea (n = 5) were calculated for the long-term average temperature of 9.7°C, following Millero and Poisson (1981). Likewise, density from salinity measurements by IOW (n = 169; 1996–2018) was also calculated (Supplementary Table 2). A density difference of 0.006 g cm⁻³ between Baltic Sea waters (IOW(Leibniz Institute for Baltic Sea Research Warnemünde), 2020) and groundwater in the Hütelmoor is more than four times lower than the value used by Robinson et al. (2007) to simulate variable-density flow experiments. As such, less convective mixing is expected due to the lower density difference of groundwater and seawater.

The effect of density on the sea level pressure head and SGD was tested by calculating the equivalent freshwater head for a water depth of 1.0 m, where most of the SGD comes out. The calculated head difference between the reference scenario and the density-corrected head is 0.01 m, with discharge flux decreasing by only 0.002 m² d⁻¹ (2%). We, therefore, assume that the model results are reliable despite neglecting density effects.

Furthermore, tides were not accounted for by the model, as they are assumed to be minimal at the study site. The Baltic Sea is an enclosed basin with a maximal tidal height of 23 cm (100-year period estimate) for the whole region and 6 cm in the nearby Danish straight (Medvedev et al., 2016).

Based on our simulations, it can be suggested that terrestrial SGD occurs more in peatlands that are 1) rewetted and diked with low peat K_s and 2) lowly-degraded peatlands with higher K_s and elevation. Future investigations on matter fluxes are thus recommended to focus on these environments.

Most coastal peatlands along the Baltic Sea are assumed to have low peat hydraulic conductivity as a result of draining and diking activities arising from agricultural use and coastal protection measures. Many of these peatlands may be rewetted in the future like the Hütelmoor since the restoration of peatlands by increasing the height of the water table is a cost-effective measure to cut down on greenhouse gas emissions (Leifeld and Menichetti, 2018). Low K_s and diking result in less infiltration of surface waters and ponding in the peatland. This water infiltrates the dune and eventually flows out as SGD. One interesting finding of this study is that most waters that come out as SGD are from the ponded surface waters.

Conversely, we can also expect SGD in peatlands that have higher elevation and are less degraded. In this setting, the peatlands have higher K_s allowing more water to infiltrate the ground, flowing into the dune and out of the seafloor as SGD. The peat extension, likely compacted due to pressure, may still hinder upflow from the aquifer sands. However, the upper peat layer will allow more groundwater to flow to the beach. Less degraded peat will also have a higher percentage of macroporosity, leading to preferential flow pathways. These observations point out the importance of land use in the past and land management activities in the future.

In both mechanisms, we have demonstrated that even a narrow dune can be important for groundwater recharge and subsequent SGD generation, as reported earlier for larger dune belts (Stieglitz, 2005; Röper et al., 2012). These results have consequences not only for the quantity but also for the quality of the discharged water, emphasizing the role of the dune and beach in total SGD flux. Taken together, our findings suggest that the dune is important for recharge and infiltration of peatland surface waters.

It is now well-known that groundwater often has higher concentrations of carbon and materials (Moore, 2010) than its receiving body of water where it discharges. Over the last few years, organic-rich subterranean estuaries have been shown to have higher concentrations of remineralized components of organic matter, including nutrients, DOC, DIC, DOM, trace elements, reduced species such as ammonium and iron (II), and lower pH (Taniguchi et al., 2019). At our study site, DOC concentrations measured in groundwater in peat and aquifer amount to 48–490 mg C L⁻¹ and 12–99 mg C L⁻¹, respectively, and to 38–380 mg C L⁻¹ in the peatland’s surface water. Our values are much higher than the previously mentioned DOC average from SGD (Szymczycha et al., 2014). While DOC may undergo several biogeochemical transformations before it is discharged to the sea, the high concentrations in peatland surface and groundwaters offer a glimpse of potential large DOC inputs via SGD. In addition, the absence of large rivers in the southern Baltic Sea could increase the importance of SGD as pathways for water and material transport.

Our study showed that a peat layer extending to the sea, and potentially cropping out, hardly influences the quantity of water fluxes. However, the peat could still be important biogeochemically. Peat from the study site has been shown to release DOC in contact with saline water (Gosch et al., 2018) and low-saline groundwater (Kreuzburg et al., 2020). Exposure of peat to water with changing salinity can promote the release of DOC and remineralized components such as CO₂ and DIC (Kreuzburg et al., 2020). As such, peat deposits along coastal zones could be potential hot spots of increased release of these materials and may be important to the release of climate-relevant gases (Kreuzburg et al., 2020). Moreover, additional geological complexities such as marine sediment bulldozing by bivalves and other groups in marine environments (Santos et al., 2012) and peat degradation and preferential solute transport on land (Liu et al., 2017) can enhance material export from coastal peatlands. Given the knowledge gaps in material transport and its abundance in the Baltic Sea, coastal peatlands warrant further scientific investigation to address their potential as an overlooked source of water and materials for the Baltic Sea via the SGD pathway, characterized by low hydraulic conductivity and low hydraulic gradient.