Drivers of seasonal- and event-scale DOC dynamics at the outlet of mountainous peatlands revealed by high-frequency monitoring

Peatlands store ∼ 20 %–30 % of the global soil organic carbon stock and are an important source of dissolved organic carbon (DOC) for inland waters. Recent improvements for in situ optical monitoring revealed that the DOC concentration in streams draining peatlands is highly variable, showing seasonal variation and short and intense DOC concentration peaks. This study aimed to statistically determine the variables driving stream DOC concentration variations at seasonal and event scales. Two mountainous peatlands (one fen and one bog) were monitored in the French Pyrenees to capture their outlet DOC concentration variability at a high-frequency rate (30 min). Abiotic variables including precipitation, stream temperature and water level, water table depth, and peat water temperature were also monitored at high frequency and used as potential predictors to explain DOC concentration variability. Results show that at both sites DOC concentration time series can be decomposed into a seasonal baseline interrupted by many short and intense peaks of higher concentrations. The DOC concentration baseline is driven, at the seasonal scale, by peat water temperature. At the event scale, DOC concentration increases are mostly driven by a rise in the water table within the peat at both sites. Univariate linear models between DOC concentration and peat water temperature or water table increases show greater efficiency at the fen site. Water recession times were derived from water level time series using master recession curve coefficients. They vary greatly between the two sites but also within one peatland site. They partly explain the differences between DOC dynamics in the studied peatlands, including peat porewater DOC concentrations and the links between stream DOC concentration and water table rise within the peatlands. This highlights that peatland complexes are composed of a mosaic of heterogeneous peat units distinctively producing or transferring DOC to streams.

1 Introduction

Aquatic carbon transfer from terrestrial ecosystems to inland waters is receiving increasing attention as it plays a major role in the watershed carbon balance (Webb et al., 2018) and in the global carbon cycle (Cole et al., 2007; Drake et al., 2017). The origin of aquatic carbon has been tracked and wetlands have been shown to be the main organic carbon suppliers to rivers at both local (Hope et al., 1997; Laudon et al., 2004; Ledesma et al., 2017) and continental scales (Hope et al., 1994; Spencer et al., 2013). Peatlands are specific wetlands which have accumulated organic matter through slow vegetation decomposition processes (Joosten and Clarke, 2002; Limpens et al., 2008). Peatlands grow under different climates (Broder et al., 2012; Dargie et al., 2017; Gorham, 1991; Page et al., 2011) and store between 20 % and 30 % of the total global soil carbon stock (Leifeld and Menichetti, 2018; Nichols and Peteet, 2019; Scharlemann et al., 2014). Stream outlets of peatlands have been monitored at different latitudes (Billett et al., 2006; Leach et al., 2016; Moore et al., 2013) in order to quantify and understand the aquatic carbon transfer between these organic-carbon-rich pools and their draining streams. Dissolved organic carbon (DOC) is a key component of these fluxes as it contributes to more than 80 % of the aquatic carbon exported from peatlands (Dinsmore et al., 2010; Hope et al., 2001; Müller et al., 2015; Roulet et al., 2007). At the outlet of peatlands, DOC is not only considered for its role in the carbon balance but also because it may be an issue for water treatment quality (Ritson, 2015) and a conveyor of potentially harmful elements along inland waters (Broder and Biester, 2017; Rothwell et al., 2007; Tipping et al., 2003).

Variability in the DOC concentration signals at the outlet of peatlands has been observed at the interannual (Fenner and Freeman, 2011; Köhler et al., 2008), the seasonal (Leach et al., 2016; Tipping et al., 2010) and even the event scales (Austnes et al., 2010; Dyson et al., 2011). DOC concentrations were found to be negatively correlated with discharge in boreal systems (Köhler et al., 2008), positively correlated with discharge in temperate areas (Clark et al., 2007) or noncorrelated with discharge in mountainous areas (Rosset et al., 2019). Temperature was also reported as an important driver of seasonal variations of DOC concentration in field (Billett et al., 2006) and mesocosm (Pastor et al., 2003) experiments since DOC production is boosted by a greater vegetation and microbial activity during warmer periods. Higher temperatures were also shown to enhance evapotranspiration from peatland resulting in a rise in DOC concentration in peat porewater and stream waters during dry summer periods (Fraser et al., 2001). Studies have highlighted that the heterogeneity of the hydraulic conductivity within peatlands (Rycroft et al., 1975) influences the water table level fluctuations (Bernard-Jannin et al., 2018; Kalbitz et al., 2002; Strack et al., 2008) and the oxygenation of the acrotelm (Freeman et al., 2001), thus driving DOC production and its transfer to streams.

DOC concentration monitoring at the outlet of peatlands has generally consisted in a weekly or monthly stream water sampling routine (Clark et al., 2008; Juutinen et al., 2013). Higher-frequency sampling has been restricted to specific high-precipitation events (Austnes, 2010; Clark et al., 2007) or snowmelt (Laudon et al., 2004). Recently, new optical in situ sensors (Rode et al., 2016) were used to track DOC concentration at a high-frequency rate (∼30 min) at the outlet of peatlands (Koehler et al., 2009; Ryder et al., 2014; Tunaley et al., 2016), highlighting the strong variability of the DOC concentration signal over a year. While diel DOC concentration cycles have been analyzed under steady hydrological conditions (Tunaley et al., 2018), no analysis has yet been performed to understand the high-frequency variability of the DOC concentration at a multiyear scale.

Mountains host many small peatland areas that are often neglected in global peatland assessments but which drastically influence stream chemistry in headwater catchments (Broder and Biester, 2015; Rosset et al., 2019). The harsh mountainous climatic conditions (from the montane to the alpine belt, Holdridge et al., 1967) and the relief of those areas generate high gradients of different abiotic parameters (temperature, precipitation, hydrology) evolving along both seasonal and event (snowmelt, rainstorms) scales. In the present study, a bog and a fen in the French Pyrenees mountains were monitored for stream DOC concentration using an optical high-frequency in situ sensor placed at their outlet. The scientific objectives of this study were (1) to statistically identify the main abiotic parameters driving stream DOC concentration variability at each site, (2) to identify the temporal scale of these drivers and (3) to compare the DOC concentration patterns of two contrasted peatlands regarding their hydrological functioning.

2 Study sites

The peatland of Bernadouze (Fig. 1b) is situated in the eastern part of the French Pyrenean mountains (42^∘48^′9^′′ N, 1^∘25^′25^′′ E). The peatland lies at 1343 m a.s.l. It belongs to a 1.4 km² watershed on limestone rocks dominated by Mont Ceint (2088 m a.s.l.) and particularly steep (average slope =50 %). From a postglacial lake, a fen developed for 10 000 years at the Bernadouze site, reaching a peat accumulation depth of 2 m on average and more than 9.5 m at extreme locations (Jalut et al., 1982; Reille, 1990). As surficial runoff contributes to the water supply of the peatland, it is considered a soligenous (minerotrophic) fen (Joosten and Clarke, 2002). The fen is subject to an oceanic climatic influence, but weather conditions can locally be contrasted due to the specific mountainous topography. For the years 2015 to 2018, the mean annual temperature was 7.9±0.3 ^∘C and the mean annual precipitation was 1797±265 mm. Subzero temperatures and snow events are regularly observed at the Bernadouze site from mid-October to mid-May with a snow cover lasting around 85 d (Gascoin et al., 2015) from December to April and sometimes exceeding 2 m in height. Beech forest is the dominant vegetation cover in the watershed, except for the highest grassland areas (>1800 m) and the 4.7 ha of the peatland (3 % of the watershed area). Vegetation on the peatland is mainly composed of species characteristic of minerotrophic peatlands such as Carex demissa and Equisetum fluviatile. However, some ombrotrophic species such as Sphagnum palustre and Sphagnum capillifolium are observed on the southern part of the peatland, forming small hummocks and revealing a progressive disconnection with the stream and the water table supply. Selective logging (one tree out of three was cut) was carried out during autumn 2016 in the lowest forested area surrounding the peatland, producing no clear hydrological and biogeochemical changes at the outlet of the peatland.

Figure 1 (a) Location map of the Ech bog (brown plot) and Bernadouze fen (green plot) in southwestern Europe. Satellite views of the peatlands of (b) Bernadouze (1343 m a.s.l.) and (c) Ech (706 m a.s.l.) and location of the site instrumentation. Map source: Esri, DigitalGlobe, Geoeye, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo and the GIS user community.

The peatland of Ech (Fig. 1c) culminates at 710 m a.s.l. in the west-central part of the French Pyrenees (43^∘4^′59^′′ N, 0^∘5^′39^′′ W). Dominated to the north by mount Cossaout (1099 m a.s.l.), the peatland depends on a 0.86 km² watershed principally composed of grasslands and grazing areas. The bog area is 5.3 ha (6 % of the watershed area) and the peat deposit reaches 3.3 m in the center (Millet et al., 2012). Peat formation started about 8200 years ago from a postglacial lake dammed by a recessional moraine in the south (Rius et al., 2012). The peatland is classified as a bog since the surface vegetation depends only on water supplied by precipitation. The site experiences a mountainous oceanic climate characterized by an average annual temperature of 11±0.2 ^∘C and an annual precipitation of 1242±386 mm (data from 2015 to 2018). Subzero daily mean temperatures are rare (∼10 d yr⁻¹) and snow events are sparse in Ech. From the model of Gascoin et al. (2015), the average duration of snow cover does not exceed 10 days at this altitude in the Pyrenees. The vegetation observed is typical of ombrotrophic bogs with a large blanket of Sphagnum capillifolium and Sphagnum compactum. Small birches and hummocks of Molinia caerulea have started to develop within the peatland. Many burning events have been reported on the peatland since its formation (Rius et al., 2012). Nowadays, agropastoral practices still use fire to limit the vegetation height and Molinia caerulea extension. The last burning event at the Ech site occurred 8 weeks before the stream monitoring in April 2017 and concerned the northeastern half of the peatland. A second burning event occurred in February 2019 in the western area of the site. It was decided to stop data acquisition just before the fire to avoid potential shifts in DOC concentration induced by this anthropogenic disturbance (Brown et al., 2015). According to field observations and a previous study related specifically to DOC exports (Rosset et al., 2019), these two mountainous peatlands are considered the main DOC contributors in their watershed.

3 Material and methods

3.1 Site instrumentation

This article presents high-frequency data monitored from 1 September 2015 to 31 December 2018 at the Bernadouze site and from 22 May 2017 to 19 February 2019 at the Ech site. Precipitation (liquid and solid) and air temperature were recorded every 30 min at Bernadouze (Gascoin and Fanise, 2018) and every 60 min at Ech by automatic weather stations located respectively 300 and 15 m from the peatlands in open areas. At both sites, sensor failures prevented data acquisition, and gap-filling models were used to complete the datasets. For missing precipitation data in Bernadouze (27 % of the monitored timeline), a linear model (r²=0.99, p value <0.01) based on cumulative precipitation recorded in Saint-Girons (414 m a.s.l.; 42^∘58^′58^′′ N, 1^∘8^′45^′′ E) was built to generate total daily precipitation. A similar model was built in Ech (r²=0.99, p value <0.01) based on data recorded in Ossen (517 m; 43^∘4^′0^′′ N, 0^∘4^′0^′′ W) to gap fill 80 % of the timeline. Missing air temperature data (5 % of the timeline) were estimated at Bernadouze from a linear regression model (r²=0.99, p value <0.01) based on data monitored at the same rate under the forest canopy 100 m away from the main weather station. In Ech, daily mean temperatures were estimated (80 % of the timeline) using a linear regression model (r²=0.88, p value <0.01) with daily mean temperature recorded in Tarbes (360 m a.s.l.; 43^∘10^′55^′′ N, 0^∘0^′2^′′ W).

At the outlet of each peatland, a multiparameter probe (Ysi EXO2, USA) measured fluorescence of the dissolved organic matter (fDOM, λ excitation $=\mathrm{365}\pm \mathrm{5}$ nm and λ emission $=\mathrm{480}\pm \mathrm{40}$ nm), turbidity, water level and temperature every 30 min. Wiper sensors prevented the optical sensors from biofouling before each measurement and the probes were inspected and calibrated monthly. In Bernadouze, battery or sensor dysfunctions and wiper failures prevented data acquisition during 14 % of the monitored period. At both sites, a network of piezometer wells (eight in Bernadouze and four in Ech) was used to record hourly the water table depth and the water temperature with automatic probes (Orpheus Mini water level logger, OTT HydroMet, Germany). Piezometer locations were selected so as to be representative of the different topographic and vegetation surfaces observed on each peatland (hummocks, lawns, river banks) (Fig. 1). The piezometer wells are 50 mm diameter PVC tubes slotted from the bottom to 10 cm below the soil surface. The average depth in Bernadouze is about 1.25±0.30 m, except for two piezometers in the center of the peatland which were drilled to 2.18±0.02 m depth. The average depth of the Ech piezometers is 2.35±0.05 m (Table S3 in the Supplement).

3.2 Water sampling and DOC calibration

Water grab sampling was performed every 2 weeks at the outlet of Bernadouze peatland and every two months at the outlet of Ech. Piezometer wells were used to sample peat porewater on four occasions (2013, 2014, 2015, 2018) in Bernadouze and on two occasions (2017, 2019) in Ech during stream baseflow periods. Water grab samples were collected using a manual peristaltic pump and were directly filtered on-site using 0.22 µm cellulose acetate filters (GSWP04700, Merck Millipore, USA). To avoid contamination from cellulose, the first milliliters of filtered water were discarded. Water samples were brought back to the laboratory in a cool box and were stored at 6 ^∘C until analysis. High-resolution water sampling was performed during 9 flood events at the outlet of Bernadouze and once at Ech using automatic water samplers (ISCO 3700, USA) to collect water during various hydrological conditions. Each sampling event consisted in collecting 24 samples of raw water (950 mL) at a frequency defined thanks to the observed time lag of discharge (1 h for rainfall and 4 h for snowmelt-driven flood events). Floodwater samples were collected within the 48 h following the previous sampling and processed as water grab samples at the laboratory.

For all samples (grab and flood samples), nonpurgeable organic carbon (NPOC, referred to hereafter as DOC) concentration was analyzed in filtered samples after acidification to pH 2 with a TOC-5000A analyzer (Shimadzu, Japan). The quantification limit was 1 mg L⁻¹. Above this value, the analytical uncertainty was estimated at ±0.1 mg L⁻¹. Reference material included ION-915 ($\left[\text{DOC}\right]=\mathrm{1.37}\pm \mathrm{0.41}$ mg C L⁻¹) and ION 96.4 ($\left[\text{DOC}\right]=\mathrm{4.64}\pm \mathrm{0.70}$ mg C L⁻¹) (Environment and Climate Change Canada, Canada).

The fluorescence of DOM (fDOM) data was explored for potential adjustments for temperature, inner filter effect and turbidity (Downing et al., 2012; de Oliveira et al., 2018; Watras et al., 2011). fDOM data were corrected for temperature as described by de Oliveira et al. (2018). The inner filter effect was adjusted at Ech for data showing absorbance values at 254 nm higher than 0.6 (de Oliveira et al., 2018). Lastly, fDOM data recorded during high-turbidity events (>20 FNU, formazin nephelometric unit at the beginning of high-discharge events were ignored in the analysis as the fluorescence can be drastically attenuated by the presence of particles (Downing et al., 2012). These periods were sporadic, accounting for only 0.2 % of the fDOM time series, and they do not alter the fDOM variability, which is delayed compared to the turbidity (Rosset et al., 2019). High-frequency DOC concentrations were calculated at each site using a site-specific linear model ($\left[\text{DOC}\right]=a\cdot \text{fDOM}+b$) linking corrected fDOM data to DOC concentration in floodwater and water grab samples. The two models (Fig. S1 in the Supplement) are respectively described by the following parameters: a=0.192, $b=-\mathrm{0.031}$, number of observations =174, r²=0.93, and p value <0.001 for Bernadouze and a=0.290, $b=-\mathrm{1.359}$, number of observations =28, r²=0.73, and p value <0.001 for Ech.

3.3 Water level fluctuation characterization

In order to provide an overall characterization of the peatlands, a mean peat water table depth, as well as a mean water temperature, was calculated at each site by averaging peat water table depths and water temperature data at a given time from the set of piezometer probes. Calculations were performed only when all sensors were running (94 % of the time period in Bernadouze and 100 % in Ech). Hereafter, the mean water temperature in the piezometers is assimilated to peat water temperature.

Master recession curve (MRC) analyses were performed on water table and stream level time series, using the MRCTools v3.1 software (Posavec et al., 2017). In order to characterize the hydrodynamic properties of the peat, MRC were preferred to hydraulic conductivity estimations from slug tests because they can be performed directly with the water table level datasets and repeated easily on other peatlands. The MRC represents the average recession of the water level observed when only discharge flow occurs (no recharge). An exponential master recession curve was used to adjust the observed average MRC and to define a specific recession coefficient (α, unit $=d^{-\mathrm{1}}$) characteristic of each monitoring point (Eq. 1) where B is a constant.

$$\begin{array}{}\text{(1)}& \mathrm{\text{Master recession curve}}\iff \text{water level}=f\left(t\right)=B\cdot e^{-\mathit{\alpha }t}\end{array}$$

The exponential recession coefficient corresponds to the inverse of the average water recession time, called recession time, in the area of a piezometer or in a stream after a precipitation event. In the following, the recession time coefficient (1∕α) is used to characterize the hydraulic properties of peatlands and stream.

3.4 DOC peak selection and characterization

Peak selections in the DOC concentration timeline were performed by running Python 3.6 (Python Software Foundation, 2019) scripts using the function find_peak available in the SciPy signal library (Jones et al., 2001) and the arithmetic mean of the DOC concentration signal (DOC_mean) as an input parameter. The peak selection criteria were as follows: to reach DOC_mean concentration and have a prominence higher than 0.25 times DOC_mean. Peaks occurring during an interval shorter than 12 h apart were grouped under the highest DOC concentration peak. Each DOC concentration peak was defined by the time period delimited by the two nearest low points surrounding the peak event. Low points were located on the DOC concentration timelines by applying the find_peak function on the negatively transformed ($-\mathrm{1}*$) DOC concentration signal previously processed with a Savitzky–Golay filter (window length =23 and polyorder =2). Low points occurring during an interval shorter than 12 h apart were grouped under the lowest DOC concentration point. Lastly, the DOC peak period could be manually adjusted to fit or correct a peculiar peak pattern. A DOC concentration peak period was characterized by different metrics (Fig. 2): its initial value corresponding to the DOC concentration of the low point at the start of the peak period (DOC_initial); its maximum value corresponding to the DOC peak value (DOC_max); and its range (DOC_increase), which was calculated by subtracting the initial value from the maximum value and finally by the rising time duration (rising_limb) which separates the initial low point time from the peaking time. In this study, initial values and increases of DOC were the targeted variables to be explained. Initial values of DOC were used to determine a DOC concentration baseline (Fig. 2). The following classification was used to describe seasonal variations: winter (December, January, February), spring (March, April, May), summer (June, July, August) and autumn (September, October, November).

Figure 2 Characterization of DOC concentration peaks. Peak events are identified on the DOC concentration timeline in blue. Each DOC concentration peak event is defined by an initial concentration (green points) and a maximum one (red points). DOC concentration increase is calculated by subtracting the initial from the maximum concentration. The time between 2 maximum DOC concentrations corresponds to the duration (seconds) separating two events and is used as an explanatory variable. The DOC concentration baseline (orange dotted line) corresponds to the time series defined by all the initial values of each DOC concentration peak.

3.5 Explanatory variable selection and characterization

In order to investigate DOC concentration variabilities (at two temporal scales: peak event and seasonal), nine explanatory variables were chosen (Table 1). Variables were calculated for each DOC concentration peak event using similar metrics to those previously described in the DOC peak characterization section (Fig. 2).

Table 1 Targeted and explanatory variable description.

The variables were abiotic parameters, chosen because they have been reported in the literature to have an explanatory potential for stream DOC concentration variability (Table 1). Two categories of variables were distinguished depending on whether the process they described was related to the production of DOC within peatlands or to the transfer of DOC from peatlands to streams. After sensitivity tests and in accordance with the observations of Tunaley et al. (2018), a mean of 7 d prior to the event was defined as the best operator to characterize pre-event conditions of air and stream water temperatures.

3.6 Correlation and statistical modeling

Relationships between targeted variables (DOC_increase and DOC_initial) and the explanatory variables were investigated using ordinary least squares (OLS) multiple regression analyses. Prior to the analyses, variables which did not satisfy a normal distribution were log or square root transformed to improve normality (Table 1). Multicollinearity was assessed among all the predictors using Pearson correlation with a threshold |r<0.7| following Dormann et al. (2013). When two variables were found to be collinear, we selected the one that displayed the highest absolute correlation with the targeted variables. Then at both sites, all variables were standardized to a mean of zero and a standard deviation of one to derive comparable estimates in the following analysis. We performed a backward stepwise selection procedure on the full model (i.e., the model including the variables retained after removing multicollinearity) to capture the best set of variables explaining each targeted variable. At each step of the procedure, the nonsignificant variables (p value >0.05) with the highest p value were dropped from the model, and the resulting reduced model was reevaluated. This process was continued until there were no nonsignificant variables remaining in the final model. To account for the time dependency of the variables in the analyses, time was also included as an explanatory variable in the full model. This variable corresponds to the duration which separates each DOC peak event from the start of the timeline. Residuals of the final models were surveyed in order to detect deviations from normality and homoscedasticity and to identify outliers. No specific deviations or outliers were detected. Model residuals were also checked for autocorrelation to verify the absence of any cyclical variation in the variables set. When more than one variable was retained in the final model, the relative contribution of each variable was assessed using hierarchical variance partitioning (Chevan and Sutherland, 1991). According to the previous predictor selections for the multiple linear regressions (MLR) models of DOC concentration increases (DOC_increase), OLS regression analyses were performed at each piezometer plot of a peatland site, replacing the mean water table increase variable by the specific water table level increase values of each plot in order to test the importance of recession time heterogeneity in the observed correlations. Similar OLS regression analyses were performed at the outlet of streams by replacing the mean water table increase variable by the stream water level increase when necessary. R² and relative importance (%) of the stream or water table level increase variable were reported for each OLS regression tested. All the analyses were undertaken in R (R Core team, 2019) using the package rms (Harell, 2019) and relaimpo (Groemping and Matthias, 2018).

4 Results

4.1 Climate, hydrology and DOC dynamics

Climatic variables are contrasted between the two studied areas. In 2018, temperatures were higher in Ech than in Bernadouze with an annual mean air temperature, water temperature and peat water temperature respectively of 11.3, 10.7 and 11.9 ^∘C compared to 7.9, 7.1 and 7.7 ^∘C. Total precipitation reached 2151 mm in Bernadouze and 1140 mm in Ech. In these steep mountainous headwaters, short and intense flood events were triggered by strong precipitation events and/or the snowmelt. Over the whole timelines, the maximum and mean of the stream water level were respectively 1.36 and 0.35 m in Ech and 0.81 and 0.10 m in Bernadouze. These short flood events were followed by recession sequences revealed by the slow decreases in the peat water table at both sites, especially in late summer and autumn (Fig. 3c). The average and minimum of the peat water table depth in the two piezometer networks were respectively −0.23 and −0.43 m at Ech and −0.15 and −0.45 m in Bernadouze. No clear relationship was observed at either site between the stream water level and the peat water table time series. The peat water levels responded differently to rain events depending on the season. For instance, a strong flood observed in the stream can be contiguous with a low or high peat water table rise (i.e., July 2016 and February 2017 events in Bernadouze) (Fig. 3b, c). Peat porewater, occasionally sampled in the piezometers, showed an average DOC concentration of 12.4±8.3 mg L⁻¹ in Bernadouze while it reached an average of 37.3±18.8 mg L⁻¹ in Ech (Table S3). Peat porewater was on average more acidic in Ech (pH $=\mathrm{5.0}\pm \mathrm{0.4}$) than in Bernadouze (pH $=\mathrm{6.2}\pm \mathrm{0.3}$).

Figure 3 Precipitation and air temperature (a), stream temperature and water level (b), high-frequency DOC concentration (c), and mean water table depth variation and peat water temperature (d). Time series observed at the outlet of the Bernadouze fen (left panels) from 1 September 2015 to 31 December 2018 and at the outlet of the Ech bog (right panels) from 22 May 2015 to 13 February 2019. The vertical gray lines represent a change of year. Green (for Bernadouze) and brown (for Ech) plots in time series (c) refer to DOC concentration measured in water grab samples and automated flood samples.

DOC concentration was highly variable at both sites during the monitored periods as highlighted by the numerous short DOC peak events (∼30 h duration) in the two time series (Fig. 3 and Table 2). At the Bernadouze site, DOC concentration peaks showing higher values were more frequent from April to November, while this was less obvious at the Ech site where DOC concentration also peaked during winter. In 2018, the arithmetic means and flow weighted averages of DOC concentration were higher at the outlet of Ech, reaching 7.1±6.1 and 4.6 mg L⁻¹, than in Bernadouze, where they were 2.0±1.5 and 1.7 mg L⁻¹.

Table 2 Time series and DOC concentration peak metrics in Bernadouze over the 1 September 2015 to 31 December 2018 period and in Ech over the 22 May 2015 to 13 February 2019 period. Mean notations correspond to arithmetic means which are given with standard deviations.

4.2 DOC concentration peaks characterization

Peak characterization (Table 2) revealed that the increases and the maxima of DOC concentration peaks were on average 2 times higher in Ech than in Bernadouze. However, the ratio between the mean increase and the mean initial value of DOC concentration was higher in Bernadouze (2.3) compared to Ech (1.9). DOC concentration peaks occurred more often at Bernadouze compared to Ech (0.24 vs. 0.16 peak per day on average), while their duration was slightly longer (32±14 vs. 28±16 h). Rising limbs of DOC concentration peaks lasted on average 10±5 and 13±14 h at Bernadouze and Ech respectively, and they were slightly longer than the stream water rising limb averages monitored at the outlet of the two peatlands. In contrast, rising limb duration of the water table in Ech was clearly longer (22±12) compared to Bernadouze (13±7 h).

General mean and seasonal means of initial DOC concentrations were 2.5 and 3.1 times higher at Ech compared to Bernadouze (Table 3). However, at both sites, initial DOC concentrations showed a clear seasonal variability. The lowest values were observed in spring and the highest in autumn, while in summer and winter DOC concentration was close to the annual mean. DOC peak event frequencies also varied at the seasonal scale (Table 3). The highest frequencies were reported in autumn at both sites. The lowest peak frequencies were observed in winter at Bernadouze and in summer at Ech.

Table 3 Reduced models explaining DOC concentration during peak events (DOC_initial and DOC_increase) at the outlet of Bernadouze and Ech peatlands. Reduced models were obtained after a backward stepwise selection procedure applied on the full model (see details in Sect. 3.6). Adjusted R² values of each model are given as the predictors and their associated coefficient, p value and R² contribution.

4.3 DOC concentration variations models

Prior to multiple regression analyses, the air temperature over 7 d, the maximum stream water level and the initial level of the water table were excluded from the analysis because of their strong correlation with other variables (Pearson's correlation |r>0.7|) (Fig. S2). Multiple linear regressions (MLRs) followed by backward stepwise selections showed that respectively 55 % and 44 % of the seasonal variation of DOC (DOC_initial) was explained by the final models at Bernadouze and Ech (Table 3). Peat water temperature was reported as an important predictor at both sites (72 % of the variance explained by the final model at Ech and 44 % at Bernadouze). In Bernadouze, variance is similarly explained by the time between two peaks (44 %). Along the 2 years of monitoring in Ech, the strong DOC concentration values observed during the dry autumn 2018 (Fig. 2) created a positive general trend in the DOC concentration baseline. This peculiar trend drastically influenced the statistical analysis, and consequently the variable time became a significant predictor at the seasonal scale. Considering the high relative importance of the peat water temperature in the two final models, two simple linear models (Fig. 4a) were built based on this variable to illustrate the seasonal DOC concentration behavior in Bernadouze (slope =0.08, intercept $=-\mathrm{0.16}$, n=231, R²=0.26, p value <0.001) and in Ech (slope =0.10 intercept =0.50, n=100, R²=0.27, p value <0.001). In the final models of increase of DOC concentration, water table increase was the most important variable at Bernadouze (67 % of the variance explained) and the single variable at Ech. In Bernadouze, other variables such as stream water temperature, stream water level increase and the time between two peaks were significant enough to be integrated in the final model of DOC concentration increase. The R² associated with the models varied strongly between the two sites, reaching 0.77 in Bernadouze and only 0.27 in Ech. Since water table increase was the main explanatory variable for the DOC concentration increase model, two simple linear models were built (Fig. 4b) with the following parameters: slope =8.44, intercept $=-\mathrm{1.06}$, n=231, R²=0.68, and p value <0.001 for Bernadouze and slope =6.39, intercept =0.84, n=100, R²=0.27, and p value <0.001 for Ech.

Figure 4 Relationships between (a) peat water temperature and natural logarithm of DOC concentration initial value and (b) square root of water table increase and natural logarithm of DOC concentration increase during peak events at Bernadouze (green) and Ech (brown). Regression coefficients (intercept and slope), p values and R² are given in each panel.

4.4 Relationships between DOC dynamics and recession time

In the fen of Bernadouze the recession times in the peat ranged from 15 to 77 d, whereas in the bog of Ech they were longer, ranging from 53 to 143 d (Fig. 5). Stream recession times were shorter at both sites, reaching 4 d in Bernadouze and 9 d in Ech. Results of the OLS regressions conducted at each peat water level monitoring plot using DOC increase final models revealed that recession time influenced the model's efficiency (Fig. 5a). Piezometers characterized by shorter recession times showed greater determination coefficients R² (Fig. 5a). Peat water table increase was the most important predictor (pie charts Fig. 5a) for all piezometer plots, contributing at least 47 % of the explained variance of the DOC increase models. In Bernadouze, the model based on stream water level was weaker (R²=0.37) than the models based on peat water table data, while in Ech the model based on stream water level was unable to explain at all the DOC increase variation (R²=0). Recession times showed a positive relationship with DOC concentration measured in the piezometers and in the streams, with higher concentrations being associated with longer recession times (Fig. 5b).

Figure 5 Relationship between water recession time coefficients and (a) the R² of the DOC_increase MLR models or (b) the DOC concentration of each water level monitoring plot at the peatland of Bernadouze (green) or Ech (brown). In both graphs, piezometer plots are represented by solid circles while the mean of the piezometers at each site is surrounded in black. Stream plots correspond to the two black striped circles in each graph. In graph (a), pie charts represent the relative importance of the water level increase variable in the R² of each model. In graph (b), a marker represents the mean DOC concentration of a plot and vertical segments the standard deviation.

5 Discussion

5.1 Long-term high-frequency monitoring

To our knowledge, this is the first time that stream DOC concentration and abiotic drivers, including peat water table depth fluctuations, have been analyzed at the outlet of peatland sites on a multiyear period at this frequency (30 min). Previously, DOC concentration variability was investigated either at lower frequencies (Clark, 2005; Dawson et al., 2011) or during shorter periods (Austnes et al., 2010; Koehler et al., 2009; Tunaley et al., 2016; Worrall et al., 2002). Recently, high-frequency monitoring of nutrient dynamics in watersheds has developed and has revealed an unexpected variability of mobilization processes for these nutrients (Blaen et al., 2017; Rode et al., 2016; Tunaley et al., 2016). These acquisitions have allowed scientists to characterize the hot moments in the biogeochemical cycles of a watershed (McClain et al., 2003). A contribution of our study is to sequence extremely brief DOC concentration peaks and to statistically disentangle their event and seasonal drivers using synchronous high-frequency monitoring of climatic and hydrological parameters. The representativeness of both seasonal- and event-scale statistical models is enhanced by the large number of events (252 peaks in Bernadouze and 101 peaks in Ech) captured at all seasons (Table 2).

5.2 Peat water temperature controls seasonal DOC concentration baseline

Clear seasonal variations in the DOC concentration baseline were observed at both sites (Fig. 3 and Table 2). The DOC concentration baseline increased in late spring, peaked in autumn, decreased during winter and reached the lowest levels in early spring. Similar seasonal DOC concentration patterns have been observed at the outlet of other peatland sites in temperate regions (Austnes, 2010; Broder and Biester, 2015; Clark et al., 2005; Tunaley et al., 2016; Worrall et al., 2006; Zheng et al., 2018) or after the snowmelt event in boreal areas (Jager et al., 2009; Köhler et al., 2008; Laudon et al., 2004; Olefeldt and Roulet, 2012; Whitfield et al., 2010).

In this study, linear regression models revealed that the seasonal variations of the DOC concentration baseline are mostly driven by peat water temperature (Table 3). At peatland sites, temperature is often identified as a DOC concentration driver at the seasonal scale (Billett et al., 2006; Clark et al., 2008; Dawson et al., 2011; Koehler et al., 2009). Warmer temperatures directly enhance DOC production by stimulating vegetation and microbial activity (Kalbitz et al., 2000; Pastor et al., 2003). Warmer temperatures are also indirectly linked to DOC production processes in temperate and northern peatlands since they often correspond to dry periods that lower water table levels. When the water table decreases, the enzymic latch (Freeman et al., 2001) is initiated on a greater volume of acrotelm (oxygenated peat) and enhances DOC production within the upper peat layers. DOC concentration relationships with peat water temperature have already been described in an acidic fen in France (Leroy et al., 2017) and in blanket peatlands from the North Pennine uplands in the UK (Clark et al., 2005); however, in these cases DOC concentrations were measured in peat porewater. A complementary study in the North Pennines (Clark et al., 2008) showed that peat porewater DOC concentrations and stream DOC concentration were strongly correlated, meaning that, by extension, the relationship between peat temperature and stream DOC concentration could be verified for these sites.

5.3 Water table increase controls DOC concentration peaks at the event scale

This study, coupling high-frequency stream DOC concentration and water table depth monitoring at both peatland sites, revealed that peat water table increase is a strong predictor of stream DOC concentration increase at the event scale (Table 3 and Fig. 4b). Until now, stream DOC concentration variability at the event scale has been investigated in terms of discharge but rarely in terms of peat water table variation. Several studies have reported stream DOC concentration increases at the outlet of peatlands during flood events (Austnes, 2010; Ryder et al., 2014; Tranvik and Jansson, 2002; Yang et al., 2015), whereas others showed dilution during high flow events (Clark et al., 2007; Grayson and Holden, 2012; Laudon et al., 2004; Worrall et al., 2002). At the outlet of peatlands, nonlinear discharge–DOC-concentration relationships have been reported (Roulet et al., 2007; Tunaley et al., 2016) and modeled (Birkel et al., 2017); this seems to be the case at our sites where stream water level explains the variability of DOC increases during flood events only poorly (Bernadouze) or not at all (Ech) (R² contribution in Table 3 and Fig. 5a). Nonlinear/hysteretic patterns (Hendrickson and Krieger, 1964; Walling and Foster, 1975) between DOC and discharge are commonly observed in upland watersheds (Jeong et al., 2012; Strohmeier et al., 2013) and are analyzed to infer DOC export mechanisms. However, these patterns cannot predict stream DOC concentration as an MLR model integrating water table increase appears to do.

The link between DOC dynamics and peat water table has been largely investigated at the seasonal scale (Kalbitz et al., 2002; Strack et al., 2008; Hribljan et al., 2014) or in mesocosm experiments (Pastor et al., 2003; Blodau et al., 2004). The peat water table is usually considered a DOC production driver as it controls the oxygenated acrotelm volume (Billett et al., 2006; Freeman et al., 2001; Ritson et al., 2017). Therefore, different studies attempted to quantify the effect of water table position on DOC production rate in peatlands. In fen and bog mesocosms, Pastor et al. (2003) observed no DOC concentration variation in the stream water after long-term peat water table decreases. Contrastingly, increasing DOC concentrations were observed during the rewetting phase of the acrotelm at fen sites in Germany (Kalbitz et al., 2002), in Canada (Strack et al., 2008) and in the USA (Hribljan et al., 2014). Clark et al. (2009) reported similar observations after rewetting peat cores in controlled laboratory conditions. Our results are in line with these studies. Moreover, thanks to the high-frequency survey, they highlight, in addition to DOC production processes, specific hydrodynamic processes driving DOC export from peatlands at the event scale.

The correlation between DOC peak and water table increase can have different origins. First, the water table increase could create a piston flow that expels prefloodwater (Małoszewski et al., 1983), enriched in DOC. At our sites, the delay (a few hours) between stream discharge peaks, peat water table increase and DOC peaks suggests that DOC concentration peaks are not directly related to water pressure. As previously observed in peat-dominated headwater catchments by Rodgers et al. (2005), this observation rejects the piston flow hypothesis. Secondly, as DOC is mostly produced in the oxygenated and unsaturated peat volume above the water table (Billett et al., 2006; Freeman et al., 2001; Ritson et al., 2017), it can be flushed by floodwater during a flood event (Boyer et al., 1997). Due to the exponential decrease in hydraulic conductivity properties (Rycroft et al., 1975), prefloodwater under the water table is less mobile than floodwater located above (Quinton et al., 2008). Our data support this second hypothesis, with a very fast increase in DOC concentrations and a rapid DOC concentration recession in the same order of magnitude as subsurface flow recession. This does not exclude the possibility that a fraction of prefloodwaters may reach the stream during the recession time, but this mixing process is minor compared to flushing processes. These mechanisms are in line with the two-layer hydrology–biogeochemistry model developed by Birkel et al. (2017) at the outlet of a peatland. Moreover, this model emphasizes the positive relationship between the stream DOC concentration and the water table connection to the upper soil horizon. Following this second hypothesis, if DOC production was the limiting factor, the linear regression (Fig. 4b) should show a plateau for the high value of water table increase. This is not the case. Thus the limiting factor appears to be the amount of water brought by floods and ultimately the full saturation of the peat. These observations support the practices for degraded peatland restoration, where a general rise in the water table is recommended to limit water table increases and the DOC concentration peaks induced at their outlets (Höll et al., 2009; Strack and Zuback, 2013).

5.4 Hydrological influences on DOC concentration at the outlet of peatlands

The higher DOC concentration observed in summer could be explained by evapotranspiration processes that concentrate solutes in stream water. However, the evapotranspiration rates in these mountainous environments are low (<300 mm yr⁻¹) compared to precipitation (>1200 mm yr⁻¹) and should not drastically influence the seasonal DOC concentration baseline. These Pyrenean peatlands are considered the main source of DOC in their watershed (Rosset et al., 2019). However, the outlet DOC concentrations are influenced by a dilution process resulting from the input of water flowing from upstream nonpeaty areas. The hydrographs recorded at the outlet of the peatlands are in favor of uniform contributions of water from the watershed to the outlet discharge. Considering this hypothesis, the different DOC concentration levels between the two sites (Table 2) may be partially explained by the different peatland coverage in the watershed of the two sites (Bernadouze 3 % vs. Ech 6 %).

In Bernadouze, DOC concentration remained extremely low when the fen was snow-covered, and it did not drop drastically during the spring snowmelt as has been observed in boreal areas (Laudon et al., 2004; Leach et al., 2016). This pattern can be explained by (1) the low initial DOC concentration, which prevents a clear dilution being observed during the snowmelt event, and (2) the snowmelt regime in this Pyrenean catchment, which may be less sudden than in boreal regions and occurs from the early snow deposit to the beginning of the growing season, continuously diluting the low winter DOC production within the peatland.

In Bernadouze, contrary to the initial hypothesis (Table 1), the time between peaks was a negative significant predictor in both seasonal and event DOC concentration models (Table 3). This is considered to be an indirect consequence of the seasonal temperature control on DOC concentration. Indeed, snow cover and the low temperatures associated with high water table positions prevent the occurrence of DOC peaks in winter, creating large time gaps between two events (Table 2) of low initial values. In contrast, DOC production is amplified during warmer periods, resulting in more frequent stream DOC concentration peaks starting at higher initial values. In Ech, where average annual temperatures are higher and snow cover is reduced, the initial hypothesis was verified since DOC concentrations were stronger in autumn after the long summer times between peaks (Table 2). However, the variable was not significant enough to be integrated in any final model (Table 3).

5.5 Contrasted DOC dynamics related to recession times

Spatial analysis of water table variation within the peatland revealed that the studied sites are composed of several peat units, characterized by contrasted recession times. In these mountainous peatlands, recession times are related to DOC dynamics, driving model efficiency between DOC concentration increase and peat water table rise and explaining DOC concentration in peat porewater.

In the present study, both stream and peat porewater DOC concentrations were higher at Ech compared to Bernadouze (Tables 2 and S3 and Fig. 5b). This is consistent with mesocosm (Pastor et al., 2003) or field (Chanton et al., 2008; Chasar et al., 2000; Moore, 1988) peat porewater observations which highlighted higher concentrations in bogs compared to fens.

Sphagnum species, which are dominant on bogs, usually produce relatively less labile and reactive DOC than vascular plants, which are more abundant on fens (Chanton et al., 2008). Lower pH values are expected to reduce DOC solubility (Clark et al., 2005). However, these relationships are not observed at our sites. As can be seen on Fig. 5b, peat porewater DOC concentrations are related to MRC, with higher concentrations being associated with longer recession times. Beyond this hydrologic control, other parameters, such as residence time and vegetation cover, linked to bog and fen conditions, influence DOC concentration levels in peat porewater.

On average the bog of Ech presented a longer recession time (111 d) than the fen of Bernadouze (20 d). However, a large variability is observed within each site. For instance, a specific unit in the fen of Bernadouze was characterized by a long recession time of 77 d. This unit shows surface bog vegetation and topographic patterns but is surrounded by typical fen units characterized by shorter recession times (Fig. 5). Thus a peatland complex must be considered a patchwork of different units and not as a uniform peat entity.

At the event scale, the univariate model between DOC concentration and peat water table increase showed a non-negligible intercept at Ech contrasting with the model of Bernadouze (Fig. 4b). This means that, in Ech, DOC concentration increases can occur without water table increases. In this case, DOC is transferred from the upper peat layers via fast runoff flows without any water table level fluctuation. Such a phenomenon is consistent with the lower hydraulic conductivities (longer recession times) measured in bogs (Fig. 5). In contrast, DOC stored in the upper peat layers of fen units is transferred to the stream by fast percolating water raising water table levels and supplying subsurface flows (Fig. 6). This explains why the DOC increase model based on peat water table increase is particularly efficient for fen units characterized by short recession times (Fig. 5a). Recession times, used as proxies of the hydraulic conductivity, also explain the differences in peat porewater DOC concentration observed between bog and fen sites. In the fen, recession times are short, meaning that the upper peat layers are rapidly washed by precipitations, inducing sudden DOC pool depletions of the peat porewater (Fig. 3c). At the bog site, DOC stored in the upper peat layers is slowly released to the stream after precipitation events and contributes to maintaining a high stream baseline (Fig. 3c) and peat porewater DOC concentrations (Fig. 5a).

Figure 6 Schematic overview of a peatland complex. Size of the arrows corresponds to DOC quantity mobilized from distinct peatland units. The DOC concentration observed in the stream depends on the contribution of the different peat units within the peatland complex.

Thus, stream DOC concentration modeling at the outlet of peatlands must account for different proportions of fen-like or bog-like units in peatland complexes to reflect the real seasonal and event DOC concentration variability. Every unit supplies DOC to the stream at a different rate depending on its volume, distance from the stream and recession time (Fig. 6). This end-member mix concurs with the model of Binet et al. (2013) describing event and seasonal water table variability in peatlands using a double porosity parametrization. In that sense, recession time appears as a physical parameter able to characterize peatland units beyond the binary typology of bog or fen. This would surely improve the efficiency of hydrological and biogeochemical models. In the case of peatland complexes characterized by long recession times, further investigations of peatland runoffs and subsurface flows are needed, analyzing denser and stream directed piezometer transects in order to build stronger DOC concentration models.

6 Conclusion

This study reports a statistical analysis of the stream DOC concentration variability at the outlet of two mountainous peatlands. Multiyear in situ high-frequency (30 min) monitoring revealed that at both sites DOC concentration time series can be decomposed into a seasonal baseline interrupted by many short, intense peaks of higher concentrations. At the seasonal scale, DOC concentration baseline variations are mainly explained by peat water temperature, which controls integrative DOC production processes within the peatland. During the hot moments of peak events, DOC concentrations are well explained at both sites by water table increases within the peatlands. Recession time is a relevant parameter to explain peat porewater DOC concentration and the different model performances observed between bog and fen sites. Recession time assessments in different locations on the two studied sites showed that peatlands are composed of different units presenting contrasted hydraulic conductivities. Thus, peatlands should not be considered to be uniform landscapes. Distinct peatland units within the same peatland complex contribute differently to the DOC transfer processes to inland waters. Recession time assessment in piezometers appears to be a simple and promising tool to investigate hydrological processes occurring in peatlands over time and space. Indeed, water table time series are often underused and only account for a seasonal mean or minimum depth. Assessing recession times on peatlands is a first step to taking peatland water table dynamics into consideration and to explaining potentially related biogeochemical processes.