Understanding natural controls on N and C biogeochemical cycles is important to estimate human impacts on these cycles. This study examined the spatiotemporal relationships between time series of weekly monitored stream and groundwater N and C (assessed by NO3 and dissolved organic C [DOC]) in the forested Wüstebach catchment (Germany). In addition to traditional correlation analysis, we applied wavelet transform coherence (WTC) analysis to study variations in the correlation and lag time between the N and C time series for different time scales. Median transit times were used to connect hydrologic and water chemistry data. We defined three stream-water groups: (i) subsurface runoff dominated locations with strong seasonal fluctuations in concentrations, short transit times, and strong negative C/N correlations with short time lags, (ii) groundwater dominated locations, with weaker seasonal fluctuations, longer transit times, and weaker C/N correlations with lags of several months, and (iii) intermediate locations, with moderate seasonal fluctuations, moderate transit times, and strong C/N correlations with short time lags. Water transit times could be identified as key drivers for the C/N relationship and we conclude that C and N transport in stream water can be explained by mixing of groundwater and subsurface runoff. Complemented by transit times and the hydrochemical time series, WTC analysis allowed us to discriminate between different water sources (groundwater vs. subsurface runoff). In conclusion, we found that in time series studies of hydrochemical data, e.g., DOC and NO3, WTC analysis can be a viable tool to identify spatiotemporally dependent relationships in catchments.

The cycles of nutrients such as N and C are closely linked. For example, N saturation is linked to C limitation in the soil for microbial processes (Kopáček et al., 2013). In aqueous ecosystems, increasing NO3 concentrations (Galloway et al., 2008) and widespread significant increases in the concentrations of dissolved organic C (DOC) in streams during the last two decades have become a major global issue (e.g., Evans et al., 2005; Galloway et al., 2008; Vitousek et al., 1997). Newly created reactive N (biologically or photochemically active) contributes to environmental problems such as air pollution, eutrophication, soil acidification, and climate change (Vitousek et al., 1997; Galloway et al., 2008). The environmental and therefore also social consequences require detailed spatiotemporal models to predict the fate of reactive N across environments (e.g., Mulholland et al., 2008; Taylor and Townsend, 2010; Barnes and Raymond, 2010). For this purpose, it might be important to understand causal links to the spatiotemporal patterns of DOC.

Dissolved organic matter consists of anthropogenic and natural organic matter. The natural component is closely linked to biogeochemical nutrient cycling, which includes different forms of N (e.g., Thomas et al., 2014; Fork and Heffernan, 2014; Taylor and Townsend, 2010). During denitrification, NO3 is reduced to N2 gas in a heterotrophic pathway. The process is coupled to the metabolic oxidation of organic matter (Canfield et al., 2010) and is a major pathway for reactive N removal and transformation into N2, which is less biologically available. In a river network, denitrification has impacts on reactive N by reducing the N load (Seitzinger et al., 2006). Taylor and Townsend (2010) found a consistent nonlinear negative relationship between DOC and NO3 across hydrological systems. Thomas et al. (2014) extended insights into the relationship between DOC and NO3 by looking at not only the quantity but also the quality of DOC. They found a negative linear relationship between DOC quality (i.e., biodegradability of DOC) and NO3 when looking at normalized ultraviolet (UV) spectra with an isosbestic point. In this study, we assumed that a non-stationary relationship exists between DOC and NO3 that varies across spatial and temporal scales. We further hypothesized that, besides biogeochemical processes, one key controlling factor is the contribution of groundwater to the stream and its tributaries, which is linked to the water transit time.

Traditionally, mathematical methods that are used for the examination of periodicities in the frequency space, such as Fourier analysis, assume that the processes have stationary behavior and therefore do not change with time (Cazelles et al., 2008; Grinsted et al., 2004). Previous observations of the dynamics of DOC and NO3 (Bernal et al., 2005; Heffernan and Cohen, 2010; Halliday et al., 2012; Lutz et al., 2012; Worrall et al., 2015), however, have shown that their temporal dynamics may vary (Taylor and Townsend, 2010; Mengistu et al., 2013; Thomas et al., 2014). Therefore, a tool is required to investigate these different non-stationary dynamics. Wavelet analysis accounts for these non-stationary relationships by detecting patterns in signals (time series) at different time scales (Grinsted et al., 2004; Cazelles et al., 2008). Wavelet analysis has been applied in fields from medicine (e.g., Keissar et al., 2009) to hydrology (e.g., Biswas and Si, 2011; Graf et al., 2014; Fang et al., 2015) and numerous other fields. The method is also considered well suited for the analysis of ecological time series (e.g., Bradshaw and Spies, 1992; Grenfell et al., 2001; Keitt and Urban, 2005; Cazelles et al., 2008; Zhang et al., 2014). Wavelet transform coherence (WTC) analysis reveals correlations between two time series that are timescale specific. In the past, there were only a few studies that included DOC and NO3 (Kang and Lin, 2007; Rusjan and Mikoš, 2010) or water quality time series in general (Mengistu et al., 2013) in wavelet analysis. Generally, a spatial component was not included. In this study, we applied WTC to extend the insights into the relationship of DOC and NO3 for longer time series data from different sampling locations within a headwater catchment.

A recent study performed in the Wüstebach catchment by Stockinger et al. (2014) gave new insights into the spatial heterogeneity of transit times for the individual tributaries that contributed to the stream, using intra-annual variations of stable isotopes δ18O (and δ2H) in sampled waters. Hence, working in this catchment allowed us to subdivide the data set into water sources with different transit times. We assumed that correlations between NO3 and DOC concentrations will differ spatially according to the water bodies to which they refer. We then tested the benefits of WTC analysis for elucidating the non-stationary coupling of DOC and NO3 cycles.

Materials and Methods

Site Description

The Wüstebach headwater catchment (50°30′ N, 6°19′ E) is located near the German–Belgian border (Fig. 1) in the German Eifel low mountain ranges and is part of the Eifel National Park. The catchment is part of the Rhenish Massif, and bedrocks in the catchment are dominated by Devonian shales with occasional sandstone inclusions (Asselberghs et al., 1936). The bedrock is generally impermeable, with a hydrologic conductivity of 9 to 10−7 m s−1 (Geological Survey North Rhine-Westphalia, 2009) but fractured and covered by a periglacial solifluctuation layer of about 1- to 2-m thickness, separated into a main layer on the top (50 cm) and a base (50–150 cm) layer, with different permeabilities for water. The hydraulic conductivity of the different soil horizons was described by Bogena et al. (2013). The base layer has a higher bulk density and thus lower hydraulic conductivity than the top layer (Borchardt, 2012). Cambisols (Inceptisols in the USDA classification) and Cambisols–Planosols have developed on the hills and hillslopes covering the main part of the catchment, whereas in and near the creek valley, more Gleysols and Histosols are found (Fig. 1). The Wüstebach catchment is a small headwater catchment of the Erkensruhr basin that flows to the north into the Obersee–Rursee, the second largest reservoir in Germany. The Wüstebach catchment is one of the research sites within the Eifel/Lower Rhine Valley observatory, which is part of the larger Terrestrial Environmental Observatories (TERENO) initiative (Zacharias et al., 2011; Bogena et al., 2012, 2015). Briefly, TERENO comprises a German network of integrated observation platforms for the long-term investigation of the consequences of global change for terrestrial ecosystems.

The investigated catchment has an area of 38.5 ha with altitudes ranging from 595 m asl in the north to 628 m asl in the south (Bogena et al., 2015; Stockinger et al., 2014). The sizes of the subcatchments that contribute to the Wüstebach main stream vary from 0.3 to 25.3 ha (Fig. 1; Stockinger et al., 2014). A smaller tributary catchment (11.4 ha), situated immediately northeast of the Wüstebach (Fig. 1), was used for comparison purposes because the catchment has similar site specifics.

The climate can be classified as warm temperate–humid, with a mean annual temperature of about 7°C (Havlik, 2002; Zacharias et al., 2011) and a mean annual precipitation of 1220 mm (period 1979–1999; Bogena et al., 2010).

Waters are low in minerals, with the specific electrical conductivity ranging between 37 and 517 μS cm−1 in the surface water and between 60 and 429 μS cm−1 in the groundwater.

The litter layer is very variable in thickness and reaches thicknesses of up to 14 cm (Borchardt, 2012). After the Second World War, the catchment was reforested, originally for timber production (Gottselig et al., 2014), with Norway spruce [Picea abies (L.) H. Karst.] and intercalations of Sitka spruce [Picea sitchensis (Bong.) Carriere] near the creek in the riparian zone (Etmann, 2009). In addition, in the wetter and partly waterlogged areas of the catchment, alder [Alnus glutinosa (L.) Gaertn.] can also be found (Deckers, 2010). Beside the coniferous forest, larger patches of open vegetation covered by herbaceous plants, fen sedges, ferns, and mosses have developed on an old earth deposit and likewise in wetter areas (Deckers, 2010). The area receives no direct influence from agriculture.

Sampling Approach

In 2009, a high-resolution spatiotemporal measurement network was established along the stream, with 16 stream water sampling points (W1–W16). One of the sampling locations, W7, was given up soon after installation because the stream in this tributary was frequently dry, leaving 15 stream sampling locations (Fig. 1). Additionally, nine groundwater sampling wells were established in the catchment (GW1–GW9). Four of these stations (GW4, 5, 6, and 7) were not considered in this study because of technical problems and incomplete data, leaving five groundwater stations (Fig. 1). The hydrochemical sampling stations were sampled on a weekly basis from 10 Aug. 2009 to 8 July 2013 (205 wk). Sampling was independent of hydrologic conditions, covering mainly base-flow conditions. Stream water sampling glass bottles were prerinsed with stream water and then filled. To get fresh and non-stagnant groundwater samples, prepumping was performed for about 3 min prior to sampling.

The samples were kept cool during transport to the research center in Jülich, where they were stored at 4°C until analysis for NO3, DOC, and other chemical parameters (see, e.g., Bogena et al., 2015; data available at http://teodoor.icg.kfa-juelich.de/overview-de). Although Stockinger et al. (2014) reported mean transit times, they actually calculated the medians of the respective transit time distributions (MedTT). We used these MedTT, which were based on isotopic measurements from 6 June 2009 to 31 Mar. 2011 of the same weekly grab samples used for the hydrochemical data. The isotopic analysis of water samples was performed using isotope-ratio mass spectrometry with high-temperature pyrolysis to analyze for δ18O and δ2H. All δ values are given against Vienna Standard Mean Ocean Water. For hydrograph simulation, the conceptual rainfall–runoff transfer function and the hydrograph separation model TRANSEP was used (Weiler et al., 2003; Stockinger et al., 2014). For further details, see Stockinger et al. (2014).

Preliminary Grouping of Sampling Locations

The stream water sampling stations were divided into three groups based on the previous analysis of their δ18O values by Stockinger et al. (2014) (Supplemental Fig. S1) and respective seasonal behavior. These types are characterized by the amount of groundwater (old water) contribution at each station:

  • Group A: W1, W2, W8, and W9 show strong seasonal variations in their stream water δ18O and are probably driven by a stronger contribution from the near-surface aquifer of the riparian zone. During times of high groundwater levels in the riparian zone (see, e.g., Bogena et al., 2015), surface runoff will also occur due to saturation overland flow and return overland flow. For simplification, this group is referred to as the subsurface runoff dominated group.

  • Group B: W4, W5, W6, W10, W11, W13, and W14 show limited seasonal variations in their stream water δ18O and are likely to be a mixture of surface and groundwater, and therefore this group is referred to as the intermediate group.

  • Group C: W3, W12, W15, and W16 show no significant seasonal variations in their stream water δ18O and are likely to be groundwater driven—the groundwater dominated group. Sampling locations W15 and W16 are located in the neighboring reference catchment.

Additional information already existed in relation to groupings of the stream waters and associated seasonal trends in DOC content (see Bol et al., 2015). However, neither DOC concentrations in the groundwater nor that of NO3 or their relationship have been studied in detail previously.

  • Group D contains the groundwater sampling locations GW1, GW3, GW8, and GW9. It is therefore referred to as groundwater.

  • Group E contains Groundwater Station GW2 with δ18O values typical of precipitation (wide range due to seasonal variations); therefore, we assume that water sampled at this location is considerably impacted by direct filtration via preferential flow.

Wavelet Analysis

We present the background of wavelet analysis followed by missing data and its handling. Explanations on how to read the wavelet plots is provided in the supplemental material.

Background

Whereas classical methods such as Fourier analysis or autoregressive models allow the identification of correlations of time series with stationary relationships, wavelet analysis overcomes this limitation by accounting for non-stationary relationships among data sets. Comprehensive background information for the applied wavelet analysis was provided by Torrence and Compo (1998) and Grinsted et al. (2004); our analysis was implemented in the package provided by Aslak Grinsted. Two classes of wavelet transforms exist: the discrete wavelet transform, which gives a compact data presentation and is mostly used for the reduction and compression of noise, and the counterpart continuous wavelet transform (CWT), which is more appropriate for the extraction of features, especially of complex time series with a lot of noise (Grinsted et al., 2004; Sinha et al., 2005), and was therefore considered suitable for this study.

Continuous wavelet transform analyzes time series in the time and frequency space and reveals the power at every point in time and for every possible period. Cross wavelet transform (XWT) analysis combines the CWTs of two time series and reveals the period with the highest joint power present at every point in time. In contrast to that, wavelet transform coherence (WTC) analysis reveals regions in this time–frequency space where both time series covary but do not necessarily have a high power. In this study, WTC was used because it accounts for (i) frequency, (ii) duration, (iii) degree, and (iv) point in time of the correlation, as well as for (v) every possible time span of the investigation period and (vi) possible lagging of two parameters analyzed simultaneously. The mother wavelet Morlet was applied because it is one of the most common continuous wavelet functions and is located well both in scales and frequency. The continuous and non-orthogonal Morlet wavelet (with dimensionless frequency ω0 = 6) is useful for feature extraction and provides a good balance between time and frequency localization. The Morlet wavelet also yields period information almost equivalent to Fourier analysis (λωt) with λωt = 1.03 s (Grinsted et al., 2004). It is described by (Torrence and Compo, 1998): 
graphic
where ψ0 is the wavelet function, η is dimensionless time, and i is the imaginary part.
The coherency of the XWT in the time–frequency space reveals areas with common power, which is not necessarily high (Grinsted et al., 2004). It is calculated according to Torrence and Webster (1999) as 
graphic
where S is a smoothing operator, s is the scale (period), Wn is the cross-wavelet transform, X is the time series x = [xn], and Y is the time series y = [yn]. The correlation is later on graphically represented by color. The coherence is the real part of the XWT, which is normalized by the power of the single time series. Therefore it can be understood as a localized squared correlation coefficient in the time–frequency domain (Grinsted et al., 2004). Like the traditional coefficient of determination, R2 accounts only for the linear part of the relationship, but under the assumption that lagging effects are eliminated. The imaginary part of the cross wavelet spectrum is used to quantify the lag and accounts for both time and frequency space. A smoothing operator (Stime) is used that has a similar footprint to the Morlet wavelet (Torrence and Webster, 1999): 
graphic
where SScale is smoothing along the wavelet scale x axis and Stime is smoothing in time.

To estimate the statistical significance level of the wavelet coherence, Monte Carlo methods were applied. In the graphical representations, 95% significance is indicated by thick black contour lines. Because the wavelet is not completely localized in time, the CWT has edge artifacts. Therefore, the cone of influence (COI) is introduced (Grinsted et al., 2004), which is shown as a bright, cone-shaped area. Broadly speaking, the COI is the area in which the results of the correlation are most reliable.

Missing Data

In general, wavelet analysis requires equidistant and continuous data sets. Water samples for our study were obtained at all stations every 7 ± 1.5 d, and the time series are therefore regarded as equidistant. However, analysis results are not available for every sampling date (e.g., insufficient discharge, climatic conditions). Missing data cannot be handled by the methodology, and row-wise removal of dates with missing results would strongly violate the abovementioned equidistance. Therefore data gaps had to be closed prior to analysis. The straight-forward method, which does not require extensive statistics, is linear interpolation. Because this is a critical procedure that can possibly create artificial significant high correlations, we excluded two stations due to >30% missing pairwise data, namely W13 and GW8, from WTC analysis. To verify the results of the remaining stations for the analysis of robustness, we used 30 sets of uncorrelated random numbers with the same mean and standard deviation as the data itself to fill the gaps. We then compared the WTC plots and also calculated the mean bias, mean absolute error (MAE) and root mean square error (RMSE) of the coherence R2 from the linear interpolated data set with those filled up with random numbers. For all stations, the bias over- or underestimating R2 was on average <0.1, the MAE was <0.25, and the RMSE was <0.3. Therefore, we are able to quantify the uncertainty of our WTC analysis, which is low. The errors increased along with an increasing number of pairwise missing data. The WTC plots with randomly filled gaps showed minor differences to the WTC plots with linearly interpolated data, and the interpretations of the plots did not change. We therefore consider the linear interpolation as appropriate for the selected data sets. Nevertheless we want to emphasize that the results of the WTC analysis are more reliable with as few data gaps as possible. To overcome this problem in the future, the package by Grinsted that we used could therefore be complemented with a modified wavelet function that allows the analysis of time series with nonsystematic gaps. In this context, approaches have been made by, e.g., Frick et al. (1997, 1998), Stoica et al. (2000), and Mondal and Percival (2010).

Method Structure

First, the dynamics of DOC and NO3 were investigated by looking at their time series individually. In addition, correlations of DOC and NO3 at observation times were calculated. Then the MedTTs were used to determine their influence on DOC and NO3 dynamics. The MedTTs represent the time it takes for 50% of precipitation water to pass the catchment. We assumed the existence of a strong relationship between the stream chemistry and the contribution of groundwater to the main stream. Variations in the DOC and NO3 relationship at various points in time in the time series are governed by changing hydrological conditions (groundwater contribution). Wavelet transform coherence was applied to gain insight into this temporal dependency of the DOC–NO3 relationship.

Results

Temporal and Spatial Changes in Dissolved Organic Carbon and Nitrate Concentrations

Time series of Group A (W1, W2, W8, and W9) showed a very distinct seasonal pattern (Fig. 2a and Supplemental Fig. S2). The highest DOC values occurred in summer, while the highest NO3 values were detected in winter. The average DOC concentrations varied between 4.79 mg L−1 (±1.59; W1) and 6.38 mg L−1 (±3.54; W8) but with outliers going up to 30 mg L−1 (outliers >2 SD, not shown). Concentrations increased downstream from W1 to W8 but with a slight decrease occurring at W9. The same pattern was observed in the average NO3 concentrations. In W1, they ranged between 0.1 and 6 mg L−1, while they were higher at W8, with the second quartile lying at about 6 mg L−1 (data not shown).

Within the groundwater dominated sampling sites (W3, W12, W13, W15, and W16) there was generally much less seasonal amplitude in DOC than in the previous group (Fig. 2c and Supplemental Fig. S3). Nevertheless, it was still present at the tributary sampling sites W3 and the reference stream site W16. The seasonal changes were more pronounced in the NO3 concentrations at all stations. Overall, W16 showed distinctly lower concentrations than the other sites within the group. The groundwater dominated sampling sites W3, W12, W13, W15, and W16 showed comparatively low values for DOC, ranging between 0.09 and 4.45 mg L−1, with W16 showing the widest range between 0.52 to 4.45 mg L−1. The NO3 concentrations were significantly higher and ranged from 1.21 mg L−1 at W16 up to 11.71 mg L−1 at W3. The intermediate type source sample locations (Fig. 2b and Supplemental Fig. S4) showed a stronger seasonality for both NO3 and DOC concentrations than the groundwater dominated group. Furthermore, the variations in DOC seemed to be in temporal concert for all the sampling stations.

When comparing the groundwater sampling locations, GW2 showed the highest DOC concentrations, followed by GW9, which had higher concentrations than GW1, GW3, and GW8 (Table 1). However, the DOC concentrations of the GW9 samples were still significantly lower than those of GW2 (Table 1). The other groundwater locations showed very low DOC concentrations of about 1 mg L−1 (Fig. 2d and Supplemental Fig. S5; Table 1). In contrast, GW2 and GW8 both showed the lowest values for NO3. Stations GW1, GW3, and GW9 displayed no clear temporal variations and constantly low DOC concentrations. Looking at the NO3 concentrations of Group D, GW1 showed the highest values, with about 11 mg L−1, while lower values were found at GW3 and GW8 at about 8 mg L−1. The values of GW9 were higher, with concentrations of about 6.5 mg L−1, but still below GW3 and GW8 (Table 1). Group E (GW2) had generally the highest DOC of up to 12.7 mg L−1, the lowest NO3 concentrations in groundwater with 5.7 ± 2 mg L−1, and showed comparatively large temporal fluctuations (Fig. 2e).

The analysis of the relationships between the DOC and NO3 concentrations (Table 1) revealed that the coefficients of determination (R2) were medium but significant (two-sided at 1%) at many stations, from 0.35 to 0.5 (Group A: W2, W8, and W9; Group B: W4, W5, W6, W10, W11, and W14; and Group C: W13 and W16). For the groundwater stations of Group D and W1, the R2 was very low from 0.08 to 0.12. On average, Pearson’s coefficients r in all groups were mostly high (negative), around −0.5 or −0.7 (Table 1).

The MedTT and average DOC or NO3 concentrations (for 4 yr) suggested a high and significant linear relationship (Fig. 3b and 3c). A negative linear relationship between MedTT and DOC was found, i.e., DOC = 6.998 − 0.00224 × MedTT (R2 = 0.8727; n = 15), whereas a positive relationship between MTT and NO3 was observed, with NO3 = 0.994 + 0.0251 × MedTT (R2 = 0.7107; n = 145). The normalized DOC/NO3 ratio was significantly (p < 0.001) different among the three identified source water groups (Fig. 3a), i.e., 2.57 ± 0.48 (n = 4) for the subsurface runoff dominated Group A, 0.61 ± 0.09 (n = 6) for the intermediate Group B, and 0.18 ± 0.10 (n = 5) for the groundwater dominated subcatchments (Group C).

The results from the correlation analysis suggests medium to weak linear anticorrelation between DOC and NO3 (Table 1). Considering local variations and lags within the time series (Fig. 2; supplemental material), it can be assumed that lagging effects in one of the parameters at various points in time lead to a lowering of the correlation coefficients. This effect could be elucidated using the WTC analysis, which eliminates possible temporal variations and lagging effects in the time series.

Wavelet Analysis

First, the terminology we use to describe the phase shifts revealed by the WTC analyses will briefly be explained, building on the schematics and explanations presented in Fig. 4 and the supplemental material. A mathematical correlation analysis allowing for a phase shift component, like Fourier cross-spectra and WTC, does not distinguish between a strongly lagged (i.e., >1/4 period) positive correlation and an anticorrelation. Because negative correlations between NO3 and DOC have been reported and discussed before see above and, e.g., Taylor and Townsend, 2010), we prefer to describe such strong lags as anticorrelation. Therefore, the lag or phase shift is always given in relation to the smallest distance to a perfectly positive or perfect negative correlation. In the following paragraphs, the type of corresponding correlation (positive or negative) is mentioned. Where a phase arrow has a left-to-right component, indicating a phase angle between −1/2π and 1/2π and thus a generally positive correlation, the exact lag in period fractions or weeks is given with respect to such a positive correlation. For other phase angles, in contrast, lags are specified with respect to a perfect anticorrelation. For example, a phase angle of −3/8 periods (Fig. 4b) at a period of 32 wk is described as “DOC is lagging 4 wk from a perfect anticorrelation” (1/8 period = 4 wk) or “Nitrate is 4 wk ahead from a perfect correlation.”

The coherence spectra of the stream and its tributaries showed a highly significant anticorrelation between DOC and NO3 (Fig. 4 and Supplemental Fig. S6). This anticorrelation occurred around the annual timescale (period lengths around 32 wk and onward), with DOC lagging between 4 and 8 wk with respect to a perfect anticorrelation to NO3. At shorter timescales, significant correlations were dominantly restricted to limited patches in time–frequency space and subject to inconsistent phase shifts. These patches either resulted from short-term processes or are artifacts of the gap interpolation in the time series. Both are possible because patches also occurred during the robustness test (see above). Therefore these short term patches are not discussed further because overinterpretation is likely.

Figure 4 shows WTC plots for characteristic stations of each group. Due to size limitations, larger WTC plots for all stations and all groups are provided in Supplemental Fig. S6.

In Group A, both W2 and W9 show disconnected patches of significant high correlations for periods <32 wk, and very high anticorrelations from periods of 24 wk (∼1/2 yr) onward, consistently until the end of the investigated periods. Nitrate is between 4 and 8 wk ahead. Phase lags in the correlation were inconsistent for the short-term correlations but mostly indicate almost perfectly negative correlations (Fig. 4 and Supplemental Fig. S6). Station W9 exhibited slightly higher R2 than observed at W2, with a time lag of 3 to 4 wk for DOC. Station W2 is very similar to W8, showing almost perfect anticorrelations (NO3 slightly ahead), and W9 is similar to W1, although at W1 NO3 is lagging by up to 15 wk and the R2 is slightly lower (Fig. 4 and Supplemental Fig. S6).

All stations in Group B had high R2 with a consistent negative lagging effect for DOC of not more than 2 to 4 wk from periods of 32 wk and onward (Fig. 4 and Supplemental Fig. S6).

The stations from Group C displayed highly variable phase shifts and patches throughout the investigated periods. All stations indicate a lagging of DOC between 4 and 20 wk toward a negative correlation, mostly from around 32 wk and onward and, for example at W3 and W12, not for the entire periods investigated. Station W16 stuck out in this group (Supplemental Fig. S6), as was demonstrated by a consistent pattern for the entire study period from 32 wk onward, comparable to W1 of Group A.

For the groundwater stations of Group D, consistent correlations were nearly absent within the COI (Fig. 4 and Supplemental Fig. S6). Station GW9, on the other hand, showed strong and consistent significant correlations for longer periods from 40 wk onward, with NO3 lagging 1/4 period. For this station, the robustness check revealed a relatively high uncertainty, meaning that correlations were significantly weaker and lower when random numbers where used to fill data gaps. The results from this station therefore are only marginally discussed. Group E (GW2) shows irregular correlation patches in the COI, but in the second half of the study period at longer periods (50–70 wk), a negative significant correlation exists with an R2 between 0.7 and 0.9 (Fig. 4 and Supplemental Fig. S6).

In general, the results from the WTC demonstrate that significant, time-lagged correlations with high coefficients of determination exist in most of the time series, which were not revealed by conventional regression and correlation analysis (Table 1). The seasonal dynamics of DOC and NO3 observed in the time series analysis are reflected by the number of weeks where the strong correlations became obvious (24–32 wk). The WTC analysis showed only small differences between the a priori defined sampling stations of Groups A (subsurface dominated) and B (intermediate). In contrast, many of the groundwater dominated (Group C) stations showed weaker anticorrelations, which, however, started at shorter periods (and sometimes ceased again at periods of approximately 1 yr). The groundwater sampling stations of Groups D and E showed the weakest anticorrelations.

Discussion

Temporal and Spatial Changes in Dissolved Organic Carbon and Nitrate Concentrations

Overall, the measured stream and groundwater NO3 concentrations in the Wüstebach catchment during the study period ranged from 0.05 to 13.5 mg L−1 in surface water and 1.37 to 13.7 mg L−1 in groundwater, with strong seasonal variations (Fig. 2; supplemental figures). This can be attributed to the forest age and the fact that catchments with spruce forests generally store and export less NO3 than those with other forest types (Emmett et al., 1993; Rothe et al., 1999; Kelly et al., 2011). Kelly et al. (2011), for example, observed from two co-located watersheds in the Fernow Experimental Forest (Norway) that a hardwood catchment exported 15 kg ha−1 yr−1, while from a Norway spruce stand nearly none was exported (0.18 kg ha−1 yr−1). This was attributed to a change in the soil substrate properties (e.g., acidification) resulting from the increasing age of the spruce vegetation. Changes in the soil pH alter the N cycle and therefore lead to a decrease in the nitrification potential of these soils.

The grouping we initially made on the basis of δ18O concentrations (values from Stockinger et al., 2014) at the surface water stations (i.e., a subsurface runoff dominated, groundwater dominated, and intermediate group; Supplemental Fig. S1), was consistent with the observed ranges and the intra-annual variations in the NO3 and DOC concentrations (Fig. 2 and Supplemental Fig. S2–S5). Local patterns in DOC and NO3 with time were in general associated with the same hydrological conditions. In the source area, all stations are either subsurface runoff dominated or of the intermediate type (up to W9) except for W3, which is groundwater dominated. This latter station is the outlet of a subsurface pipe from a former groundwater storage reservoir (Fig. 1). Water from this pipe forms the W3 tributary of the Wüstebach, which is therefore groundwater dominated, mixed with subsurface runoff (Bol et al. [2015] and field observations). All intermediate stations are located in the stream, with W5 being the only exception as a tributary. Therefore it is the only tributary affected by both ground and subsurface runoff. This might be due to specific local heterogeneities at this location. The lower NO3 and higher DOC content at W16 compared with W15 is surprising because W15 is the upstream counterpart of W16 in the reference stream. Therefore, either some additional dilution or denitrification must take place to account for the decrease in NO3. The opposite effect can be assumed for the DOC concentration, with a concentrating effect of enrichment from additional sources rather than only W15 as a source (Newcomer et al., 2012; Barnes et al., 2012; Arango et al., 2007).

Furthermore, we could also verify the previously determined two subgroups in the groundwater samples based on the signature of the DOC and NO3 trends. One group is likely to result from leakage water (GW1, GW3, GW8, and GW9: Group D) and the other (GW2; Group E) from bank filtration or direct infiltration due to preferential flow paths or a very high groundwater table up to ground level (Fleckenstein et al., 2010; Boano et al., 2014; Ward, 2016). Wiekenkamp et al. (2016) showed that preferential flow processes are common in the Wüstebach catchment. Surprisingly, though, another groundwater well (GW9), which was located on the opposite side of the stream but also in the vicinity of the stream, does not show the same DOC or NO3 seasonal features as GW2, despite being affected by similar conditions. This is corroborated by the fact that the δ18O signatures of GW2 and GW9 are also dissimilar, and the δ18O range of GW9 is more similar to other groundwater wells in the leakage water group (Group E; M. Stockinger, personal observation, 2015).

Both Wüstebach tributaries in the lower part of the catchment (below 600 m), namely W12 and W13, are groundwater dominated (Fig. 1; Table 1). Also, the stations of the comparison catchment (W15 and W16) are dominated by groundwater. An explanation for this is the shallow groundwater in the study area, which could be observed during weekly groundwater sampling and other field work (e.g., soil sampling). Observations of the water table revealed that the groundwater table is shallow in the Wüstebach catchment but closest to the surface in the lower catchment part (field observations). This probably explains why the upper catchment is only minorly affected by groundwater, while the lower Wüstebach and its tributaries are generally groundwater dominated. The results derived from the current study concerning the DOC and NO3 characteristics of the respective catchment parts and sampling locations (Table 1; Fig. 2 and 3 and Supplemental Fig. S1–S4) support these observations in the catchments’ hydrology of groundwater behavior and influence on water chemistry.

In conclusion, the upper catchment is less affected by groundwater. Station W3 is an exception, as for most of the study period it was an artificial outlet of a water reservoir fed by groundwater (see Bol et al., 2015).

Overall long-term trends in DOC and NO3 were not observed, either in surface water or in groundwater (Fig. 2 and Supplemental Fig. S2–S5). Bol et al. (2015) observed a 4-yr lowering of specific UV absorbance in headwater tributaries W1 and W8, which means a reduction in aromaticity of the dissolved organic matter and increased biodegradability (Evans et al., 2005; Rowe et al., 2014). Therefore, any indications of changes in DOC or NO3 concentrations with the regard to global change or any other process could not be observed due to the comparatively short time span of the observations. However, the TERENO platforms are intended to monitor the systems (arable cropping, grassland, or forest) for at least 20 yr as the impact of climate change may take time to materialize in key environmental parameters. Nevertheless, strong seasonal variations could be determined, especially in the surface water data. The DOC concentrations were higher in summer than in winter, with NO3 behaving opposite. Similar observations, with a negative relationship between DOC and NO3, have been made in other studies, i.e., Piatek et al. (2009), Taylor and Townsend (2010), and Williams et al. (2011).

For the groundwater dominated waters (Group C), the transit time was in general longer than for the subsurface runoff dominated stations (Group A). Figures 3a and 3b support the hypothesis of our study that generally the DOC, NO3, and DOC/NO3 ratio in Wüstebach stream waters is mainly controlled by hydrological mixing processes, while in the soil it is additionally controlled by biogeochemical cycles. These findings indicate that the DOC/NO3 ratio works as a robust indicator for the water pathways in catchments that control the composition of the water leaving the catchment.

The results from the relationships between MedTT and DOC, NO3, and DOC/NO3 (Fig. 3a and 3b) indicate that the longer the median transit time for each single raindrop, the lower the DOC concentration and the DOC/NO3 ratio and the higher NO3 concentrations will be in the stream water samples. The impact of hydrology and therefore transport pathways on DOC (Jansen et al., 2014) and NO3 concentrations have been shown in previous studies (e.g., Oeurng et al., 2010; Zarnetske et al., 2011a, 2011b; Van Gaelen et al., 2014). These are very closely linked to biogeochemical processes, which can act longer on the water the longer it remains in the soil. The MedTT is just an indicator of how long the water needs to pass the catchment. It does not provide information about what happens to the water while being inside the catchment. According to Kirchner (2016), seasonal tracer cycles can give more accurate information about young water fractions in streamflow (up to 3 mo of age) than transit times and should be considered in future studies. They could also be included in the wavelet analysis for a comparison with DOC and NO3 concentrations and their ratio at every point in time for different time scales.

Wavelet Analysis

The wavelet coherence analysis accounted for lagging effects as well as different strengths and signs of correlation at different timescales. It revealed high correlations at seasonal to annual timescales for all subsurface runoff dominated and intermediate stations (Fig. 4 and Supplemental Fig. S6). Other hydrology-related studies using wavelet coherence also found the most consistent correlations near the annual timescale (Lauzon et al., 2004; Graf et al., 2014). In both studies, correlations at shorter timescales were intermittent and generally even degraded at the shortest timescales. In our study, the latter feature was less expressed, possibly as a result of the weekly measurement frequency (Stockinger et al., 2016). The studies by Lauzon et al. (2004) and Graf et al. (2014) were based on daily averages of automatically measured, non-chemical parameters, and Graf et al. (2014) found a notable improvement in the global correlation between water budget terms when aggregating to a weekly timescale.

In our study, we identified that the first type of station (Group A, subsurface runoff dominated) showed predominantly (R2 = ∼1) negative correlations (Fig. 4 and Supplemental Fig. S6). Group B showed patterns comparable to those of Group A. Groundwater and groundwater dominated stations (Groups C, D, and E) showed significantly lower and less consistent correlations. A recently published study by Thomas et al. (2014) addressed the hypothesis of Taylor and Townsend (2010) that stoichiometric requirements during microbial joint turnover of N and C sources are responsible for a frequent, global anticorrelation between NO3 and DOC. Their results coincide with the results from the WTC analysis conducted in this study, although using a different approach. They found that the behavior of the normalized UV spectra is able to reveal an otherwise unseen isosbestic point (HIP) at 225 nm in various samples. The HIP is the point where the UV spectra of a sample set cross at a certain wavelength (Thomas et al., 2014). It was then possible to link the HIP to a simple linear relationship between DOC and NO3 using a linear model. The proportions were found to depend on the sampling location as well as the environmental conditions. The simple linear model provided by Thomas et al. (2014) is not very feasible for time series data because it has to be applied to each sample separately. The wavelet coherence analysis, on the other hand, provides the possibility of finding hidden linear correlations in a large number of data points at the same time, even if they vary in time. Nevertheless, the two studies show that the entirely mathematical approach of the WTC conducted in this study would benefit from further insights into the chemistry of DOC and NO3; in turn, the HIP method could benefit from wavelet analysis when looking at a larger time span of observations.

Therefore, this study nicely agrees with the observations made by Taylor and Townsend (2010) that there is clearly an anticorrelation between DOC and NO3, which they found across a wide range of environments (tropical, temperate, boreal, and arctic regions). A negative though nonlinear relationship between DOC and NO3 was also observed in numerous other studies, e.g., by Thingstad et al. (2008), Zarnetske et al. (2011a, 2011b), Trimmer et al. (2012), Wickland et al. (2012), Sandford et al. (2013), Thomas et al. (2014), and Ali et al. (2015). Compared with the traditional correlation analysis shown in Table 1 and discussed above, WTC analysis revealed a stronger negative correlation if only large (seasonal and larger) timescales were accounted for and lags of several weeks allowed for. This was particularly true for Groups A and B (strongly subsurface runoff affected waters). From these discrepancies between global and WTC-based correlation analyses, we can hypothesize on several constraints on the processes leading to anticorrelation between NO3 and DOC, such as the stoichiometric requirements of microbial decomposition of both (Taylor and Townsend, 2010). First, because the anticorrelation was not perfectly in phase in many parts of the space–time–frequency continuum of our study, future research should examine whether hypotheses on such microbial processes are consistent with the possible existence of time lags between increasing (or decreasing) NO3 and decreasing (or increasing) DOC. Second, the vadose zone where subsurface runoff is happening appears to be a hot spot for such processes. At sites dominated by local leakage, at least in our study, negative correlations may be more confounded by processes leading to positive correlation, such as, e.g., simultaneous dilution or uptake from sources of both components. Finally, and most prominently, negative correlations appear to dominate at seasonal and larger timescales, as opposed to week-to-week variations. This could be due either to the time constant of reaction of microbial turnover to resource availability or, more likely, to a dominance of the abovementioned processes causing positive correlation at shorter timescales.

For example, Baurès et al. (2013) demonstrated a relationship between DOC and NO3 that was positive, and thus opposite compared with our study, during high flow rates. Therefore, further studies should consider not only base flow conditions but also event water. To apply the WTC analysis time series to events, sampling at a high frequency (e.g., 5 min) would be required. As stated by Kirchner et al. (2004), weekly or monthly water quality sampling would be “like trying to understand a Beethoven symphony if one could only hear one note every minute or two!”

Conclusions

Our study revealed that DOC, the DOC/NO3 ratio, and NO3 concentrations in Wüstebach stream waters can be explained by hydrological mixing processes. The DOC concentration and the DOC/NO3 ratio increased while the NO3 concentrations decreased along with the MedTT of the stream water in the catchment. The contribution of groundwater (≈MedTT) dominated the spatiotemporal dependent relationship between DOC and NO3. The results from the wavelet transform coherence analysis reflect the different water characteristics, expressed in the DOC and NO3 relationships. Strong negative correlations between DOC and NO3 for subsurface runoff dominated sampling sites and those of the intermediate type sampling sites were revealed (R2 = 0.8–1.0). Good negative correlations were also found for groundwater dominated sampling sites (R2 = 0.7 and 0.9) when the time series were shifted by 1 to 5 mo. For groundwater sampling sites, correlations were significantly weaker and inconsistent.

The WCT analysis revealed strong negative correlations between DOC and NO3 that were consistent across longer timescales, while they were weaker and inconsistent at shorter timescales. This fact might help to realign previous findings of other researchers on different correlations at varying scales, e.g., Taylor and Townsend (2010) on global negative correlations and Baurès et al. (2013) on a positive correlation at short timescales during high flow rates.

Future studies should focus on covering all temporal scales from 5 min to years to gain insights into the timing of hydrological and biogeochemical processes.

Acknowledgments

We thank TERENO (Terrestrial Environmental Observatories), funded by the Helmholtz-Gemeinschaft. We also thank Werner Küpper, Rainer Harms, Willi Benders, Leander Fürst, and Ferdinand Engels for the weekly hydrochemical sampling and Holger Wissel for stable isotope analyses of water samples.

Open access.
This is an open access article distributed under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).