Patterns and dynamics of dissolved organic carbon exports from a riparian zone of a temperate, forested catchment

Export of dissolved organic carbon (DOC) from riparian zones (RZs) is an important, but poorly understood component of temperate catchment carbon budgets. This paper delineates explicit DOC source zones within the RZ of a small forested catchment in central Germany, and identifies and quantifies their dominant DOC export mechanism at high spatiotemporal resolution. Stream water DOC samples from differing hydrological situations were compared to riparian DOC 15 groundwater and surface water samples and classified chemically (via Fourier-transform ion cyclotron resonance mass spectrometry) and spatially via a small-scale topographic analysis of the RZ at a resolution of 1m. Explicit water fluxes from the resulting riparian DOC source zones were then simulated by a physically-based, fully-integrated numerical flow model (HydroGeoSphere). Chemical classification revealed two distinct DOC pools (DOCI and DOCII) in the RZ. The comparison of stream and riparian 20 water samples indicated a predominant export of DOCI during wet conditions and high groundwater levels. The two DOC pools were spatially separated and mapped using a threshold value in high-resolution topographical wetness index (TWIHR). Hydrological modelling revealed that surface runoff from DOCI source zones with high TWIHR values dominated overall discharge generation and therefore DOC export. Although corresponding to only 15 % of the area in the studied RZ, the high TWIHR zones provided in total 1.5 times the load of DOC from the remaining 85 % of the area associated with the DOCII pool. 25 Our results suggest that surface DOC export can play a dominant role for DOC export in RZs with overall low topographic relief and should be considered in DOC export models. We propose that proxies of spatial heterogeneity (here: TWIHR) can delineate the most active riparian source zones and provide a meaningful basis for improved model conceptualization of surficial DOC export.

An improved understanding of the dominant mechanisms in small-scale landscape elements could help to find accessible proxies that can describe the larger-scale effective DOC-export behavior of catchments (Frei et al., 2012;Grabs et al., 2012).
Proxies based on landscape-scale characteristics like different land use types (Pisani et al., 2020), hydromapping based on 70 convergence of topography Ploum et al., 2020), or topographical wetness (e.g. represented by the topographic wetness index TWI (Beven and Kirkby, 1979)) can potentially provide information on the distances and dominant connections of DOC source locations to the stream and thus DOC export at catchment-scale Fellman et al., 2017;Andersson and Nyberg, 2009;Inamdar and Mitchell, 2006;Grabs et al., 2012).
However, landscape-and catchment-scale proxies cannot differentiate DOC export mechanisms that are induced by the small-75 scale heterogeneity of topography and hydrological properties (e.g. micro-topography driven variable source zone contributions or spatial variability in biogeochemical DOC mobilization) because they integrate the entire RZ into a few spatial entities (e.g. model cells). For example, Andersson and Nyberg (2009) proposed a skewed linear relationship between mean catchment-scale TWI and DOC concentration for various Swedish catchments, which showed poor performance during dry hydrological states. We argue that a smaller-scale, dynamic assessment of the TWI, which includes the wetness conditions 80 in the RZ could improve the quality of the DOC-TWI relationship by Anderson and Nyberg, as it would be able to account for the actual contributing source zones of DOC. We postulate that a small-scale assessment of the TWI in the RZ can identify dominant DOC source zones and their export dynamics, and in turn will yield better correlations between TWI values and instream DOC concentrations than catchment-scale TWI assessments.
Several powerful tools to assess topographical, biogeochemical and hydrological heterogeneities in riparian zones have 85 become available over the years. High-resolution drone-based digital elevation models (DEMs) of study sites are now relatively easy to obtain and provide a more detailed perspective on small-scale topographical variations in the RZ that can be linked to potential DOC source areas and source apportionment. State-of-the-art techniques for molecular DOC characterization based on FT-ICR-MS allow to compare chemical fingerprints of DOC from different riparian source areas with the integrated DOC signal in the stream (Raeke et al., 2017;Seifert et al., 2016;Wagner et al., 2019).  properties of the DOC compounds can provide valuable information on the origin, state of processing and mobilization potential of DOC, which can complement the identification of riparian DOC export patterns derived from topographic analysis. Last but not least, fully integrated, hydrogeological modeling can provide detailed insights into spatio-temporal patterns of small-scale hydrologic processes in the RZ (Frei et al., 2010). A systematic combination of these methods will ultimately lead to an improved understanding of the dominant DOC source areas in RZs and the mobilization and hydrologic 95 export of DOC from these areas to the stream. https://doi.org/10.5194/hess-2021-82 Preprint. Discussion started: 23 February 2021 c Author(s) 2021. CC BY 4.0 License.
In this paper we identify dominant source areas of DOC within a riparian zone based on topographic analysis and DOC fingerprinting and evaluate the temporally variable export of DOC from these source areas to the stream using numerical flow modeling. To this end, we mapped and classified riparian DOC source areas and assessed DOC export to the stream within a highly instrumented RZ. More specifically, (1) chemical DOC fingerprints from different locations within a riparian zone were 100 compared to the integrated DOC fingerprints in stream runoff under base flow and event flow conditions. (2) The DOC source areas identified in (1) were then mapped in space using high-resolution TWI to derive hot spots of DOC mobilization within the RZ. Finally, (3) we quantify the hydrological contribution of the delineated DOC source zones from (2) to stream flow generation and solute export using a numerical surface-subsurface water flow model. This approach allows us to develop a robust proxy to explicitly identify probable DOC source areas in riparian zones under different hydrological conditions, which 105 could be used for mechanistically sound conceptualizations of catchment-scale, parsimonious DOC export models.

Rappbode Catchment
Measurements were conducted in a headwater catchment of the Rappbode stream (51°39'22.61"N 10°41'53.98"E, Fig. 1) located in the Harz Mountains, Central Germany. After draining into a drinking water reservoir, the Rappbode stream flows 110 into the river Bode, and discharges (through the rivers Saale and the Elbe) into the North Sea. The catchment has an area of 2.58 km² and a drainage density of 2.91 km km -1 . The study site is characterized by a temperate climate (Kottek et al., 2006), with a long-term mean air temperature of 6.0 °C and mean annual precipitation of 831 mm (Stiege weather station 12 km away from the study site, data provided by the German Weather Service DWD). The uncultivated and uninhabited catchment is predominantly forested with spruce and pine trees (77 %), 11 % is covered with grass, and 12 % is covered by other vegetation 115 and a few unpaved roads. Elevation ranges from 540 to 620 m above sea level; the mean topographic slope is 3.9°. The 90th percentile of the topographic wetness index (TWI) as a measure for the extent of riparian wetlands in the catchment  is 10.10 (median 7.58). The geology at this site consists mainly of graywacke, clay schist and diabase (Wollschläger et al., 2016). Soils in the spring area are dominated by peat and peat formation. Overall, 25 % of the catchment soils are humic and stagnic gleysols that are connected to riparian zones. 120 https://doi.org/10.5194/hess-2021-82 Preprint. Discussion started: 23 February 2021 c Author(s) 2021. CC BY 4.0 License.

The monitoring program
The study site was chosen to be within the 90 th percentile of the Rappbode catchment TWI (derived from a drone-based high resolution DEM) and is thus regarded to be representative for riparian sites of the catchment . Electric resistivity tomography (Resecs DC resistivity meter system, Kiel, Germany) was applied at two transects using a Wenner alpha configuration with an electrode distance of 0.5 m in order to explore structural consistencies of the subsurface. The intensive 130 monitoring campaign took place from 12 April 2017 until 19 December 2018, the overall monitoring period lasted until 23 July 2019.

In-stream data sensors
Two PCM4 portable flow meters (Nivus, Germany) measured discharge in the Rappbode stream at a chosen inlet (PCM4in) and outlet of the study site (PCM4out), respectively. A pressure transducer (Solinst Levellogger, Canada) was installed in the 135 center of the study site. All three probes measured water level every 15 minutes. Discharge at the center was then estimated via a stage-discharge relationship (stage was measured using a pressure transducer) which was established based on biweekly https://doi.org/10.5194/hess-2021-82 Preprint. Discussion started: 23 February 2021 c Author(s) 2021. CC BY 4.0 License. manual discharge measurements using an electromagnetic flow meter (n = 37; MF pro, Ott, Germany). Maximum discharge measured manually was 0.22 m³ s -1 on 28 February 2017 at a water level of 37.9 cm. Manual measurements were recorded between 28 February 2017 and 19 December 2018. The extrapolation of the stage discharge relationship to a wider range of 140 stages was found to be in a valid range (Werner et al., 2019). However, values larger than 0.22 m³ s -1 are more uncertain than smaller values.

Weather station
A weather station (WS-GP1, Delta-T, United Kingdom) was placed about 250 m northwest of the study site in order to characterize ambient weather conditions. Air temperature, humidity, wind direction and speed, solar radiation and rainfall were 145 recorded at a 30 min interval. Potential evapotranspiration (ETP) was calculated from the weather data after Penman-Monteith (Allen et al., 1998) also at a 30 minute resolution.

Piezometer network
The elevation map derived from the drone flight ( Figure 1b) revealed that the floodplain's slope in the direction of the stream (0.2m/10m) was steeper as the slope towards the stream (0.1m/10m). We expected that this would have ramifications for the 150 direction of the slope of the groundwater level and its temporal dynamics. To have maximum ability in capturing the magnitude and direction of this slope as a function of time a piezometer network was installed aligned on a square grid, with one principal axis oriented in parallel to the stream and the other perpendicular to the stream. Adjacent to the Rappbode stream, 25 partly screened piezometers (2.54 cm diameter, HDPE, 10 cm screen) were installed in a rectangular grid pattern, comprising a piezometer network covering 3600 m² (60 m x 60 m, Figure 1b). The well spacing was 12.5 m in both principal directions of 155 the grid. In addition 3 more wells were installed at 0.3 m depth inside the rectangular grid for surface near sampling. The A horizon in the piezometer holes was 17.7 cm ± 2.4 cm on average (n = 27) ( Figure 1). Figure S1 gives the depth of slotting and the soil horizon accessed by the slotted section. The screen depth of the piezometers ranged between 20 cm and 107 cm below ground (average = 75.2 cm). This was a tradeoff between having continuous water level measurements from the pressure transducers and covering the anticipated large variety of different DOC types in different soil layers and depths (Shen et al., 160 2015). Each piezometer was equipped with a pressure transducer (Levellogger, Solinst, Canada and Diver, van Essen, Netherlands), measuring at a 15-minute interval. All pressure transducers were barometrically corrected and adjusted to manual measurements of the groundwater level at 8 occasions during the 15 months measurement period (from 04 October 2017 to 19 December 2018).

Sampling and maintenance 165
Biweekly routine samples were collected in the Rappbode stream to determine DOC concentration. Riparian wells (n = 28), respective stream water samples and -if possible -riparian surface water samples were taken on five occasions throughout the year (Table S1). In addition, during five discharge-generating events, sampling in the stream was conducted hourly to sub-https://doi.org/10.5194/hess-2021-82 Preprint. Discussion started: 23 February 2021 c Author(s) 2021. CC BY 4.0 License. hourly by auto-samplers (6712 Full-Size Portable Sampler, Teledyne ISCO, US), which were triggered by the rate of water level increase. Bottles from the auto-sampler were soaked for 48 h in 0.1 N HCl prior to use. Process blanks with deionized 170 water were prepared to correct for eventual contaminations during field work and sample processing. Due to the remoteness of the study site, auto-sampled stream water samples were collected within 4 days after the triggered event sampling. Samples were stored in the dark inside the sampler and air temperatures were below 10°C during that time.
Riparian zone shallow groundwater samples were collected from 3 to 18 out of the 28 installed piezometers depending on hydrological conditions. Before sampling, water in the wells was replaced one to three times (based on the responsivity of the 175 wells) through pumping. The flasks and the pump were rinsed with sample water prior to actual sampling. 100 mL of sample volume was then transferred into acid-rinsed (0.1 N HCl) and baked (500 °C, 4 h) glass bottles and stored dark and cool until further processing in the laboratory.

Sample processing and DOC determination 180
After sample collection at the study site, samples were filtered using 0.45 µm membrane filters (cellulose acetate filter, Th. Geyer, Germany) and acidified to pH 2 (HCl, 30 %, Merk, Germany) for subsequent DOC measurement and extraction. Filters were rinsed with 20 mL of sample water to avoid bleeding. Filtered samples were stored in the dark at 4 °C until laboratory analysis was conducted (typically within two days).
DOC concentration was determined as non-purgeable organic carbon with a high-temperature catalytic oxidation system (multi 185 N/C 3100, Analytik Jena, Jena, Germany) from acidified samples and extracts after solvent evaporation. Due to the small difference in DOC concentration between hourly and subhourly samples during events we chose a lower time resolution for high resolution mass spectrometry measurements. A volume of 10 -200 mL (n = 142) was extracted via solid-phase extraction using an automated system (FreeStyle, LC Tech, Obertaufkirchen, Germany) on 50 mg styrene-divinyl-polymer type sorbens (Bond Elut PPL, Agilent Technologies, Santa Clara, CA, United States) to desalt the sample for subsequent DI-ESI-MS 190 according to Dittmar et al. (2008) and (Raeke et al., 2017). The carbon-to-sorbens ratio (C:PPL) was 280 ± 130 (m/m, n = 142). The SPE-DOM was eluted with 1 mL methanol (Biosolve, Valkenswaard, The Netherlands), and stored at −20 °C until measurement. Carbon based extraction efficiency was (56 ± 15) % (determined from n = 133 samples, Fig. S8). This is in the range of typical extraction efficiencies obtained for freshwater samples (Raeke et al., 2017). Immediately prior FT-ICR-MS analysis extracts were diluted to 20 ppm and mixed 1:1 (v/v) with ultrapure water Merck,Darmstadt,195 Germany).
France) instrument was used in ESI negative mode (capillary voltage: 4.2 kV) using an Apollo II source. Extracts were 200 analyzed in random order with an auto sampler (infusion rate: 10 µL min −1 ). For each spectrum, 256 scans were co-added in the mass range 150 -1000 m/z with 4MW time domain (resolution @ 400 m/z was ca. 483000). Mass spectra were internally re-calibrated with a list of peaks (247 -643 m/z, n > 55) commonly present in terrestrial DOM and the mass accuracy after linear calibration was better than 0.13 ppm (n = 142). Peaks were considered if the signal-to-noise (S/N) ratio was greater than four. Raw spectra were processed with Compass DataAnalysis 4.4 (Bruker Daltonics Inc., Billerica, MA, USA). SRFA 205 reference sample and a pool sample (mix of randomly picked DOM extracts) was repeatedly measured to check instrument performance across multiple measurement days and solvent and extraction blanks were measured with the samples.

FT-ICR-MS data processing
Molecular formulas were assigned to peaks in the range 0-750 m/z allowing for elemental compositions C1-60H0-122N0-2O0-40S0-1 with an error range of ± 0.5 ppm according to Herzsprung et al. (2020). Briefly, the following rules were applied: 0.3 ≤ H/C 210 Koch et al. (2014)), -10 ≤ DBE-O ≤ 10 (Herzsprung et al., 2014), and element probability rules proposed by Kind and Fiehn (2007). Isotopologue formulas ( 13 C, 34 S) were used for quality control but removed from the final data set as they represent duplicate chemical information. All peaks present in the instrument blank and in the SPE blanks were subtracted from the mass list. Relative peak intensities (RI) were calculated based on the summed intensities of all assigned peaks in each sample. To ensure that the 215 variance in DOC quality observed by FT-ICR-MS was not induced by systematic instrumental shifts at different times of the year, we quantified the variability of peak intensities based on 15 reference samples (SRFA) of the four measurement days.
Subsequently this variability was applied to every RI in measured samples to derive a mean error for the intensity (see S1).
We conclude that the analytical uncertainty from the FT-ICR-MS measurements between the different measurement dates does only minor affect the overall variance of the samples, which allows the joint evaluation of all samples (see S1, Fig ). An assigned molecular formula is termed compound throughout this article although they potentially represents multiple isomers.

Numerical water flow modeling
The numerical code HydroGeoSphere (HGS) was used to quantify flows at the study site. HydroGeoSphere is a 3D numerical model describing fully coupled surface-subsurface, variably saturated flow (Therrien et al., 2010). It solves Richards' equation 225 for 3D variably saturated water flow in the subsurface domain, and uses Manning's equation and the diffusive-wave approximation of the St. Venant equations to simulate surface flow in the 2D surface domain and 1D channel network (Yang et al., 2015). Using a dual node coupling approach, HydroGeoSphere simulates the water exchange fluxes between the domains, providing the simulated infiltration/exfiltration fluxes. More details on the governing equations, coupling approach, and general aspects of HydroGeoSphere can be found in (Therrien et al., 2010). 230 https://doi.org/10.5194/hess-2021-82 Preprint. Discussion started: 23 February 2021 c Author(s) 2021. CC BY 4.0 License.
Only the upper 2 m of the alluvial sediments were included in the flow simulation as an aquifer because geological survey data showed that the electric resistivity dropped sharply below that depth indicating base rock formation (Fig. S4). The subsurface was discretized into 8 horizontal element-layers, each composed of 6924 prisms. The layer thicknesses ranged from 0.05 m near the land surface to 0.5 m near the aquifer bottom. The horizontal cell sizes varied from 1 m to 2 m. The 6924 uppermost 2D triangles of the 3D prismatic mesh were used to discretize the surface domain. The channel crossing the study site was 235 discretized into 148 1D segments, which coincide with the segments of the 2D triangular mesh. The line element made up of the channel segments was treated as a Cauchy boundary with the stream stage being calculated based on the assumption of a rectangular cross-section with channel width and depth based on measurements.

Parameters
Lateral variability of the saturated hydraulic conductivity K, i.e. the K field in xy-plane, was calibrated using 38 pilot-points 240 (Tang et al., 2017;Moeck et al., 2015) distributed inside the study site. Each of these pilot-points were associated with a K value, set to 0.1 m d -1 prior to calibration. For the vertical K heterogeneity, it was assumed that the K was depth-dependent and decreased exponentially when the aquifer was deeper than 0.2 m, as K = K0 when d < 0.2 m, and K = K0e -λd d > 0.2 m, where K0 is the hydraulic conductivity of the aquifer top determined from the horizontal K field, is the aquifer depth to land surface and λ is a factor constraining the decreasing rate, set to 0 prior to calibration. These formulations captured the general 245 decreasing trend of K with depth, while also reflecting the fact that this decreasing trend was not significant in the upper 0.2 m of the soil, which contained most roots and mainly consisted of poorly decayed organic materials.
The surface domain and channel domain were uniformly parameterized with Manning roughness coefficients (Manning et al., 1890), respectively. Prior to calibration, the roughness coefficients were set to 6•10 -6 d m -1/3 , a typical value for floodplains/grassland. The parameters described above were selected as key parameters that could significantly influence the 250 flow processes, and were optimized during calibration (Table S2) (Table S2).

Boundary and initial conditions
Input data was defined at a one-hour time resolution for the simulation, and all the time-resolution data (15 min) was aggregated 255 accordingly. For the aquifer top boundary, spatially uniform and temporally variable precipitation was applied to the surface domain. Spatially uniform and temporally variable potential ET, estimated using the climate data, was specified as model input with actual ET being simulated by the model (Therrien et al., 2010). For the upstream boundary AB (Figure 2) were estimated using • , where is the annual mean groundwater recharge rate in this area (~200 mm yr -1 ), and is the contributing surface area associated with each lateral boundary estimated from the DEM. The respective recharge fluxes Qleft, and Qright were calculated as 0.18 m 3 s -1 per unit length and 0.09 m 3 s -1 per unit length. They were also allowed to vary by 265 0.1 to 10 times of their initial values during model calibration (Table S2)

Calibration 280
Transient calibration was performed using the software package PEST, which uses the Marquardt method to minimize a target function (describing the error between modeled and measured variables) by varying the values of a given set of parameters until the optimization criterion is reached (Doherty and Hunt, 2010). The calibrated model parameters (Table S2)  The time-variable groundwater levels were well replicated by the model for the wells near the channel (Fig. S2a). The wells close to channel had better fits than those near the side boundaries (Fig. S2b), because the latter were more strongly constrained by the constant groundwater fluxes through the side boundaries. The calibrated flow model was used to quantify internal water flux data from specific regions in the riparian zone. Additionally, advective-dispersive particle-tracking was used on the flow 295 field from the calibrated model to visualize surface and subsurface flow paths through the model domain. The surface flow paths are used to identify key runoff generation zones in the riparian zone.

Statistical methods
Statistical analysis was performed using R (R-Core-Team, 2017). Evaluation of geospatial properties was conducted via R in combination with ArcMap (ESRI, US). 300

Chemical classification of potential DOC source zones
Peak intensity weighted average (wa) of FT-ICR-MS derived molecular parameters (mass (mz), elemental ratios (H/C, O/C, N/C, S/C), nominal oxidation state of carbon (NOSC) and aromaticity index (AI)) was calculated for each sample by Eq. (1): where wap(x) is the weighted average value for the molecular parameter p in sample x. pi(x) is the derived value for the 305 parameter p of each molecular formula i in sample x. Accordingly, inti(x) is the peak intensity for molecular formula i in sample

x.
A principal components analysis (PCA) was then performed with the RI of molecular formulas (n = 482) commonly detected in all riparian samples (n = 66) covering on average 40 % of the assigned intensities in each sample. A consecutive k-means clustering on the first two principal components (R package FactoMineR (Lê et al., 2008)) was used to partition the riparian 310 https://doi.org/10.5194/hess-2021-82 Preprint. Discussion started: 23 February 2021 c Author(s) 2021. CC BY 4.0 License. samples into two (as suggested by the silhouette index (Rousseeuw, 1987)) chemically distinct groups, representing different DOC quality in the riparian groundwater.
The Wilcoxon rank sum and the Kolmogorov Smirnov (KS) test were applied to identify significant differences in the distributions and medians of DOC concentration and FT-ICR-MS derived molecular parameters for the two groundwater DOC groups and stream water samples. 315

Hydromorphological classification of potential DOC source zones
Pearson correlation analysis was applied to every groundwater level time series with the stream water level time series.
Geomorphological analysis was conducted via the TWI, according to Eq. (2) log , Where TWI is the topographical wetness index for each cell, f is the flow accumulation (the accumulated weight of all cells 320 flowing into a downslope cell at the surface) at each cell and s is the slope in radians of respective triangular surface element.
The DInf algorithm was used for calculating flow accumulation in ArcMap since it proved to depict more realistic hydrological routing (Tarboton, 1997). To account for mathematical infinity/indefinite terms, zero slopes were set to 0.001 rad, and cells with no flow accumulation (f = 0) were set to 1 cell instead.
A smoothed map of the local TWIHR values was created by assigning the median TWI value of the central cell and its 8 325 surrounding cells to the central cell. According to KS and F-test statistics, the resulting map represented the non-smoothened TWI distribution of the study site (pKS = 0.33; pF = 0.76). We applied the Wilcoxon rank sum to test for differences in TWIHR distributions and medians of the two DOC clusters. For an extrapolation from point sources (i.e. sampled piezometers) to spatial entities, zones were demarcated that had higher TWIHR values than the median of the DOC cluster group of higher TWIHR. The water balance for the entire model site and the two TWI-generated zones was then estimated and compared to 330 each other between 12 April 2017 and 19 December 2018 by modeling with HydroGeoSphere.

Hydroclimatic conditions and DOC chemical characterization
The basic statistics of discharge, groundwater level and climatic variables throughout the 15-months measurement period are given in Table 1. Discharge shows event-type, erratic variability but in general followed a clear seasonal pattern, with lowest 335 values in late summer and highest values in spring (Figure 3a). Stream water level was highest during a flood event from 01 to 03 January 2018 when the Rappbode stream went over-bank. We decided to not include this event in the statistics, because we could not estimate the discharge for water levels higher than the stream banks. Yet observing this flood event helped to verify and understand riparian surface runoff pathways at our study site. The amount of precipitation during 2018 (580 mm) https://doi.org/10.5194/hess-2021-82 Preprint. Discussion started: 23 February 2021 c Author(s) 2021. CC BY 4.0 License. was below the long-term annual mean (831 mm) at the nearest official weather station. Air temperature exhibited seasonal 340 patterns and was above the long term annual mean at the nearest station (8.6°C vs. 6.0°C).  An overview of hydroclimatic data for the dates of DOC sampling (Figure 3b, c) in the stream or the RZ is given in Table S3. DOC concentrations in the stream during events generally followed the hydrograph, with higher concentrations during higher water levels. Although DOC concentration and molecular properties across all riparian groundwater samples exhibited variability (Table 1)

Chemical classification
A detailed overview of the FT-ICR MS results of the distinct water samples can be found in S2, Table S5, and Figs. S9-S10.
The PCA to classify riparian DOC quality was able to explain 66.3 % of the total variance of DOC molecular peak intensities 375 using 2 principal components (PCs). K-means clustering based on the PCs then separated the riparian samples into two groups of 19 and 47 samples (DOCI and DOCII, resp.; Fig. S6), representing distinct DOC quality in the riparian zone. Weighted average molecular parameters were significantly different between the two clusters, allowing for a clear separation between DOCI and DOCII (Figure 4a). Samples clustered in DOCI had higher DOC concentration and their molecular composition was characterized by more oxidized (higher NOSC and waOC), more aromatic molecules (higher waAI), with a lower fraction of 380 heteroatoms (smaller waSC, waNC not shown), and a lower molecular weight (smaller wamz). Comparison of DOCI and DOCII molecular composition and concentration with that of stream water sampled during rain events in spring, summer and autumn confirmed overall different DOC quality distributions and medians between riparian groundwater and stream water (Fig. 4).
However, median values of the DOCI cluster were always closer to the median of stream water event samples than the respective DOCII median. Moreover, the DOC composition of one event in December (Fig. 4a, orange dots) was in the range 385 of the riparian samples, but did not show much compositional variability within the event.
DOCI samples from April (n = 9, Fig. 3) and December (n = 9) did not show significant differences in DOC molecular composition (except waHC) and concentration (Figure 4b). In addition, DOC concentration and quality in the stream samples (from the routine measurement program, non-event conditions) generally matched DOCI concentration and quality in April and December (except waHC and waAI). In contrast, DOCII samples from April (n = 13) and December (n = 33) differed 390 significantly according to their wamz, waHC, waAI, waNOSC and DOC concentration (Figure 4c). While DOC concentration and quality of stream water samples from December were mostly within the range of the respective DOCII samples, stream water samples from April were mostly outside the range of the DOC properties and concentrations of the respective DOCII samples.
The cluster DOCI was associated with groundwater sampled at depth to water table (DTW) > -0.3 m in 9 cases whereas the cluster DOCII had 7 samples no deeper than 0.3 m. Median high-resolution TWI (TWIHR) values at the well position (see 2.5) 395 were grouped according to their attribution to the DOCI and DOCII clusters based on the chemical characterization. Note that 8 wells, sampled during different occasions throughout the year occur in both DOC clusters and according TWIHR values can thus occur in both clusters (Fig. 4d). In general, the median values of TWIHR for DOCI wells were significantly higher (Wilcoxon rank sum p < 0.008) than respective values of the DOCII wells ( Figure 4d). Both, distribution of TWIHR was different and median was higher when comparing DOCI vs DOCII TWIHR values from April, whereas December samples did not show 400 any statistical significant difference in their TWIHR distribution or median. https://doi.org/10.5194/hess-2021-82 Preprint. Discussion started: 23 February 2021 c Author(s) 2021. CC BY 4.0 License. Table 1

Spatial mapping
The significant difference in TWIHR median values of DOCI and DOCII wells (Wilcoxon rank sum p < 0.008) was used to 415 spatially separate both potential source zones from each other by using the median TWIHR value of the DOCI group (9.66) as a threshold. Overall, this TWI threshold also includes 25 % of the samples from the DOCII group. However, this percentage dropped to 15 % in April. Also note that different samples of one well can appear in both DOC groups. This TWIHR classification split the riparian zone into zones of high TWIHR (DOCI source zone) and low TWIHR (DOCII source zone) values ( Figure 5). The high TWIHR zones defined in this way made up 14.6 % of the study site. The HydroGeoSphere (HGS) model 420 was then used to quantify the runoff generation from the delineated DOCI source zones and to quantify their impact on total runoff generation and DOC export from the study site. According to our simulations surficial runoff (that is groundwater discharging to the surface or direct precipitation onto saturated areas feeding the stream) was the main contribution to overall runoff generation at the site ( Figure 6). The median contribution of surficial runoff to total runoff generation was 61 % (± 12 % standard deviation) but surface contributions increased up to 99 % during event situations. We selected the subsurface-425 surface exchange flux as a key descriptive variable for potential surface runoff contributions, because it quantifies the availability of water at the surface for each cell of the model. Although there was a 1.5 times higher net surface water flux generation from low TWIHR zones throughout the modeling period (Figure 6b), the median of the area-normalized water exchange flux for high TWIHR zones (0.026 m d -1 ) was about 8.6 times higher than that for DOCII source zones (0.003 m d -1 ).
This resulted in higher absolute exchange fluxes in high TWIHR zones in about 47 % of the modeling period. During (non-430 winter) runoff events, water exchange flux contribution of high TWIHR zones increased up to 100 % (negative or no exchange flux for DOCII source zones in dry summer) whereas low TWIHR zones contributed more potential surface runoff at non-event winter conditions and flooding events when high overall exchange fluxes occurred under fully saturated soil conditions. Hydrological conditions were exemplary mapped for situations on 13 December 2017, immediately after an event under wet conditions, and on 29 August 2018, amidst a prolonged dry period in summer (Fig. S7, Table S4 for according water fluxes). 435 High TWIHR zones had the highest exchange fluxes and water depths in both wet and dry situations. Surface flow paths in winter intersect the high TWIHR zones establishing hydrologic connectivity between these zones and the stream, which is in line with our observations on DOC quality. The highest positive exchange flux values (GW exfiltration) occurred at the outer hillslope boundaries of the RZ, running parallel to the channel (values at the exact boundaries of the model have not been taken into account due to potential boundary effects). These exfiltration spots were located close to the strongest surface water 440 infiltration spots. During the exemplary wet situation, the entire RZ was saturated with water besides the stream banks. Surficial runoff pathways then connect the DOCI source areas (high TWIHR zones), running parallel to the stream and eventually entering it within the modeled domain.

Surface DOC export from high TWIHR and low TWIHR source zones
During the model period, DOCI source wells had a median DOC concentration of 5.8 mg L -1 which was 2.3 times higher than for the DOCII source wells. We assumed the DOC concentrations to stay in a range of mean ± SD throughout the year (cf. 455 Figure 4b, c). DOC export was then roughly calculated by multiplying mean ± SD of DOCI and DOCII concentrations with the absolute surface runoff volumes from the respective high and low TWIHR zones. With that mean overall export from high TWIHR zones exceeded that from low TWIHR zones in about 70 % of the time although making up only 14.6 % of the total area. In absolute numbers, high TWIHR zones exported roughly 1.5 times the amount of DOC (7.1ꞏ10 6 g) to the stream than low TWIHR zones (4.6ꞏ10 6 g). This amounts to a nearly 20 times higher area-normalized DOC export from high TWIHR zones 460 than from low TWIHR zones. Highest disparity between the export of the two source zones was during events in autumn and spring when water in the low TWIHR zone infiltrated rather into the ground (no DOC export from low TWIHR zones, Figure 7) while high TWIHR zones exported DOC (positive spikes). Infiltrating conditions for the high TWIHR zone only occurred during summer events when DOC export was generally at the minimum (mean daily export rates of 3.1 and 17.3 g d -1 for low and high TWIHR zones, respectively) whereas equally high DOC export occurred in winter (234.2 g d -1 and 267.2 g d -1 for low and 465 high TWIHR zones, respectively). The median export from the high TWIHR zone was above that from low TWIHR zone in nonwinter conditions, with the highest disparity between medians in spring and autumn (Figure 7). Then high TWIHR zones exhibited exfiltrating conditions, whereas water in low TWIHR zones kept infiltrating.

Chemical and hydrological classification of riparian DOC source zones
We chemically classified riparian groundwater samples with regards to DOC composition throughout the year and compared 475 them to the DOC composition in event-and base-flow stream water in order to detect dominant riparian DOC source zones during different hydrological conditions. Clustering the riparian samples revealed two distinct DOC pools in the RZ (DOCI and DOCII), differing in their concentrations and molecular characteristics, thus representing different biogeochemical and physicochemical settings within the RZ. The DOCI pool reflects processed plant-derived organic matter with low bioavailability (higher degree of aromatic, oxidized compounds, low fraction of heteroatoms). The molecular character of the 480 DOCI indicates low organic matter turnover presumably due to oxygen limitation even in the topmost soil layers, similar to typical wetland sites (Tfaily et al., 2018). Accordingly, DOCI source zones are connected to constantly low DTW values and high TWIHR values indicating waterlogging. In turn this may explain the conservation of DOCI due to anaerobic conditions (LaCroix et al., 2019). In contrast, the DOCII quality can be attributed to increased microbial processing of organic matter from organic rich top-soil layers. Matching our observations, the DOCII compounds are generally characterized by microbial, 485 secondary metabolites (aliphatic, heteroatom enriched molecules) and DOC which is not adsorbed to mineral phases (small, low aromatic and oxygen depleted molecules) as typically found in deeper soil layers (Shen et al., 2015;Kaiser and Kalbitz, 2012;LaCroix et al., 2019).
The close agreement between DOCI characteristics and according stream water samples in April and December indicates a predominant connectivity of the DOCI pool with the stream during wet conditions and high groundwater levels. Here the DOCI 490 quality did not show significant differences between April and December indicating a replete DOC pool with constant contribution to the overall DOC quality in the stream. In contrast, the DOCII composition was reflected in the stream water composition in December but not in April, indicating the influence of seasonality on this pool. Less organic matter input and lower biogeochemical process rates in winter and at the same time increased DOC export (at higher groundwater levels) may specifically deplete the DOCII pool (Werner et al., 2019). The composition of DOCII samples in April thus may represent a 495 pool of DOC of low solubility with low sorption affinity which was not depleted during high groundwater levels in winter (saturated and with larger molecular weight).
Variations in stream DOC composition also appeared at the shorter event time scale. Fully saturated riparian conditions caused hydrological mixing of the DOCI and DOCII pool in December, leading to a stream DOC signal that was in-between the two riparian DOC pools. On the other hand, DOCI and stream water DOC compositions are converging during events at non-500 saturated soil conditions. Consequently, changing source contributions from both riparian DOC pools are induced by the prevailing hydrological situation and therefore have to be considered as a general mechanism of DOC export in the catchment.
Observed shifts between distributions of stream and riparian DOC composition might be due to near-and instream processing (Dawson et al., 2001;Battin et al., 2003), but also inter-and intra-annual variability of hydroclimatic drivers like seasonality or antecedent soil conditions (Werner et al., 2019;Köhler et al., 2009;Strohmeier et al., 2013;Futter and de Wit, 2008 we showed direct links between six major DOC molecular properties and the DOC concentration of riparian and stream water samples. In comparison to other studies, which use integrated or indirect signals to derive information on DOC characteristics, like (specific) UV absorption and spectrophotometric slopes values (Ledesma et al., 2018b;Werner et al., 2019), or conductivity and pH (Ploum et al., 2020), this allows a spatial explicit alignment of riparian DOC sources zones to stream water samples at high credibility. 510

Hydrologic controls and quantification of riparian DOC export
Additionally, the physically-based model HGS independently shows that the dominant runoff generating mechanism (and thus potential DOC export pathway) at the study site is surface runoff from DOCI source zones, because respective high TWIHR zones have a different hydrological setting that produces more runoff per unit area and reacts more direct to precipitation than DOCII source zones (with low TWIHR values). A first estimation reveals an overall dominance of DOCI export from high 515 TWIHR zones during events, despite making up only about 15 % of the total study site. Consequently, DOC can be regarded as a site specific multi-facetted tracer that can provide a deeper insight into hydrologically controlled contributions of variable source zones: The whole RZ contributes rather uniformly to surficial DOC export during fully saturated soil conditions (e.g. during winter -as long as there is enough DOC available for transport) suggesting equal contributions from DOCI and DOCII source zones. Here a median large scale TWI will be enough to describe DOC export. In contrast, there is no surficial DOC 520 export during dry situations. However in-between those two extremes, event-induced surficial DOC export is a function of soil wetness which regulates biogeochemical DOC processing and potential hydrological surface export (Werner et al., 2019).
During such intermediate event situations, surface DOC export from high TWIHR zones increases whereas low TWIHR zones still depict zero surface export and stream water will shift to a DOCI dominated composition.

Implications for DOC export modeling 525
Recent studies conclude that lateral DOC export is not well researched, but is an important component of the global DOC budget (Zarnetske et al., 2018;Wen et al., 2020). In this regard, we found that surficial DOC export dominated overall lateral DOC export in our study site. Yet surface DOC export is underrepresented in current model-conceptualizations of lateral DOC export (Dick et al., 2015;Ledesma et al., 2018a;Bracken et al., 2013;Ploum et al., 2020). Reasons for the exclusion might be due to the high complexity of representing the spatio-temporal heterogeneity of surface export in modeling concepts or just 530 because it does not play an important role in other catchments. With this study on TWIHR, we present a proxy of small scale heterogeneities in surface runoff generation and respective DOC export that bears the potential to improve existing DOC export models and could lead to new approaches and concepts for better DOC export modeling. Adding a TWIHR based dynamic source zone activation term could greatly improve the mechanistic basis of lumped, parsimonious models potentially leading to a more accurate upscaling of DOC export from RZs, especially for RZs with minimal slopes that tend to have highest 535 groundwater levels and thus surface runoff. A general, coarse-scale relationship between soil moisture and potential surface runoff generation has already been proposed based on a catchment-scale topography-driven runoff proxy (Gao et al., 2019, https://doi.org/10.5194/hess-2021-82 Preprint. Discussion started: 23 February 2021 c Author(s) 2021. CC BY 4.0 License. Birkel et al., 2020). The mechanistic connection between TWI and surface DOC export in our study represents a similar general mechanism that is applicable for the whole riparian zone in small catchments (similar TWIHR values should result in similar runoff generation). In contrast, small scale topographical proxies of DOC export based on the presented TWIHR have the 540 potential to more accurately represent distinct zones of DOC export during differing hydrological situations. Then weather data can subsequently be used to estimate which distinct source zones (hot spots) contribute to DOC export during a precipitation event (hot moment), thus overcoming the restrictive dichotomy of the hot-spot hot-moment concept (Bernhardt et al., 2017).

Conclusions 545
Chemical classification of riparian groundwater samples via ultra-high resolution FT-ICR-MS revealed two distinct DOC pools (DOCI and DOCII) in the riparian zone. Degrading plant material presumably contributes most to an aromatic, oxygenrich DOC pool with high concentrations, located in regions of high wetness and local topographic depressions. This DOCI pool tends to be available for microbial degradation (Mostovaya et al., 2017), can be photodegraded (Wilske et al., 2020), and can be relatively easily removed through sedimentation (Dadi et al., 2017) or during drinking water treatment (Raeke et al., 550 2017). However, it also has a potential for disinfection byproduct formation (Wang et al., 2017) when not removed sufficiently.
The second pool (DOCII) reflects microbially processed, mobile DOC with lower concentration and larger compositional variability across seasons. Respective source zones of DOCI and DOCII can be separated and mapped by a threshold value in high-resolution TWI (TWIHR). The identification of source zones was achieved via independent measures (unsupervised chemical classification, and TWI-based physical flux modeling) indicating high credibility. Additionally, hydrological 555 modeling revealed that the dominant runoff generation mechanism in the study site was surface runoff. Here DOCI source zones, which make up 15 % of the total study site provided 1.5 times more DOC for export than the remaining 85 % of the area associated with the DOCII pool. Furthermore, highest discrepancy between the DOCI and DOCII surface export was during events at intermediate wetness states (neither completely saturated nor very dry). Overall, this is a strong indication that DOCI sources rather than DOCII sources get exported into the stream during event situations. We therefore conclude that certain 560 thresholds in TWIHR, which are based on actual wetness state, can identify explicit source zones of surficial DOC export from riparian zones. Thus TWIHR in turn defines the relative contributions from different source zones with more unique DOC source signals in the stream during dry and transitional periods whereas mixed signals occur during very wet conditions. In contrast to other studies in northern till catchments (e.g. Ledesma et al. (2018b)), this study highlights that surface DOC export from the riparian zone plays an important role for lateral DOC export from hydromorphic soils with overall low topographic relief. 565 Therefore we want to emphasize that surface export should be acknowledged in respective DOC export models. Delineating activated source zones for DOC export by topographic proxies of lateral, spatial heterogeneity (here represented by TWIHR) can help to identify source zones in existing DOC models or provide a mechanistic basis for improved model conceptualization for lateral DOC export modeling. https://doi.org/10.5194/hess-2021-82 Preprint. Discussion started: 23 February 2021 c Author(s) 2021. CC BY 4.0 License. Data availability. All datasets used in this synthesis are publicly available via the following link: https://doi.org/10.4211/hs.b32ba184414e475ba36a0bb193866ef1.
Supplement. The supplement related to this article is available online at: Author contributions. BJW, JHF, OJL, AM and GHdR planned and designed the research. BJW performed the statistical analysis and wrote the paper with contributions from all co-authors. JY edited the hydrological modelling. UW implemented 575 geophysical investigations. RG provided post-processed drone altimetry data.
Competing interests. Gerrit de Rooij is a member of the editorial board of the journal.