Evaluation of Legacy Forest Harvesting Impacts on Dominant Stream Water Sources and Implications for Water Quality Using End Member Mixing Analysis

Forests are critical water supply regions that are increasingly threatened by natural and anthropogenic disturbance. Evaluation of runoff-generating processes within harvested and undisturbed headwater catchments provides insight into disturbance impacts on water quality and drinking water treatability. In this study, an extensive hydrologic dataset collected at the experimental Turkey Lakes Watershed (TLW) located on the Canadian Shield was used to quantify sources of stormflow in legacy clear-cut (24-years post harvesting) and forested (control) headwater catchments using an end member mixing analysis (EMMA) model. Stream water, groundwater, soil water, and throughfall water quality were evaluated during spring snowmelt, stormflow, and fall wet-up. Groundwater chemistry was similar to stream water chemistry in both catchments, suggesting that groundwater is a major contributor to stream flow. The water chemistry in small wetlands within the study catchments was comparable to stream water chemistry, suggesting that wetlands are also important contributors to stream flow. Differences in wetland position between the legacy clear-cut and control catchments appeared to have a greater influence on source contributions than legacy harvesting. Results from this study provide insight into runoff-generation processes that reflect event/seasonal flow dynamics and the impacts on water quality.


Introduction
Forested landscapes are critical for the provision of safe and reliable drinking water to a majority of Canadians and 60% of the U.S. population, and 60% of the world's largest cities also rely on forests for water [1][2][3][4]. The ecosystem services of tropical and temperate forests have been valued at USD 23.32 trillion per year [5,6], with savings to global drinking water treatment costs valued at USD~4.1 trillion per year [7]. Landscape disturbance threatens water supply through alterations to hydrologic and biogeochemical cycles that ultimately change water quantity and quality; these threats are further exacerbated by compound disturbance (e.g., heavy precipitation following wildfire) [1,8,9]. Forest disturbance typically intensifies the water cycle, resulting in more water moving through the landscape [10,11], which increases the transport of sediment and nutrients to stream channels [9,[12][13][14]. Sediment and nutrients are transported downstream [15][16][17] and challenge drinking water treatment operations [1]. However, the varied hydrologic responses of forested catchments to landscape disturbance [18][19][20][21] make predicting disturbance effects difficult [22]. Explaining this variability requires an understanding of the impacts of landscape disturbance on water movement and storage in a variety of forested landscapes [20,22,23] and is essential to evaluate potential risks to drinking water treatment operations [1,22].
Harvesting of trees is a common anthropogenic disturbance in many forested ecosystems [24]. In addition to industrial forestry and silviculture, harvesting is increasingly advocated as a techno-ecological nature-based solution (TE-NBS) [25] for mitigating climatechange-exacerbated landscape disturbance (e.g., wildfire) threats to drinking water treatability [25,26]. Critically, however, harvesting can impact water storage and flow paths. Vegetation removal associated with harvesting can reduce rates of precipitation interception and transpiration and increase the water content of forest soils [11,27,28]. Additionally, heavy equipment can change the physical characteristics of forest soils by increasing soil bulk density and reducing soil porosity and infiltration rates [29][30][31][32]. These changes may increase soil saturation [33,34] and reduce the water-storage capacity of soils after harvesting [35,36]. Reduced water-storage capacity can increase the proportion of water transported via overland flow [29] and shallow subsurface pathways [37,38] that result in more frequent, quick flow events [36,39]. Notably, these combined changes result in higher annual and seasonal water yields [10,40] and increased peak flows [11].
Long-term experimental watersheds such as the Turkey Lakes Watershed (TLW) in Ontario, Canada, provide detailed historical hydro-chemical data that can be used to evaluate complex hydrological processes [41]. Previous work at TLW has shown that shallow soil water is a major contributor to stream flow during snowmelt [42] and that perched water tables over less permeable/impermeable soil layers direct flow through the upper soil profile [43]. Threshold relationships were observed between quick flow and precipitation [36], and stream flow generation is controlled through transmissivity feedback (rapid stream flow increases as the water table approaches soil surfaces) [38]. A forest harvesting experiment conducted at TLW showed that harvesting increased the proportion of water moving through the upper soil profile [10,38], reduced mean water travel times [23,37], and decreased the root-zone-storage capacity, resulting in a reduced precipitation threshold required for quick flow production [36]. Accordingly, increased instances of annual and seasonal water yields draining harvested catchments were observed, with some signs of recovery 15 years after harvesting [10]. Collectively, this work described key impacts of forest harvesting on hydrologic response during the first 15 years post harvest. However, few studies have been conducted to evaluate the longer-term (>20 years) legacy impacts of forest harvesting on hydrologic processes that control stream flow generation. [10,13,44].
End member mixing analysis (EMMA) is a commonly used tool that attempts to explain observed stream water chemistry composition as a mixture of different source areas or "end members", assuming a linear mixing process [45][46][47][48]. End members are pre-defined hydrologic source areas such as groundwater, soil water, precipitation, and snowmelt that are hypothesized to contribute water to stream flow. The dimensionality of stream water chemistry datasets is removed using principal component analysis (PCA), and potential end members are projected into the mixing subspace defined by the stream water chemistry [46,49]. End members are then evaluated for their ability to explain stream water chemistry in the mixing subspace, and mass balance equations are used to calculate relative end member contributions to stream flow [46,48]. Diagnostic tools using eigenvector and residual analysis have been developed to estimate the number of end members required to explain the stream water chemistry without prior knowledge of the existing end members [50].
EMMA has provided numerous insights regarding runoff-generation processes and water quality. Many studies define end members differently, making direct comparisons difficult [51]. However, other studies have demonstrated that three end members are often sufficient to explain stream water chemistry [48,52] and that the prominent end member contributing to stream water can change depending on antecedent soil moisture conditions [48,53]. For example, at the Panola Mountain Research Watershed in Georgia, USA, end member contributions to stream flow changed from outcrop runoff to riparian groundwater with increasing antecedent moisture conditions [53]. Similarly, the number of end members contributing to stream flow increased from two to three under higher antecedent moisture conditions, as the transient saturated groundwater end member was only contributing to stream flow under wet antecedent moisture conditions in the Kiryu Experimental Watershed, Japan [52].
Additional insights into how changing end members contributes to water quality have also been assembled. For example, snowmelt was reported to drive episodic acidification within stream water by diluting base cations at the Hubbard Brook Experimental Forest in New Hampshire, USA [54]. This dilution effect shifted during summer storms as the end member contributions from upper soil horizons became more prominent and supplied organic acids to the stream channel [54]. At the H.J. Andrews Experimental Forest in Oregon, USA, EMMA was used to identify lateral subsurface flow originating from the organic and upper soil horizons as the dominant end member contributing dissolved organic carbon (DOC) and dissolved organic nitrogen (DON) to the stream network during the rising limb of the hydrograph [55]. These studies highlight the utility of EMMA for discerning end member contributions to stream water under different antecedent soil moisture conditions and subsequent impacts on water quality. However, few studies have utilized EMMA to evaluate landscape disturbance impacts on runoff-generation processes or inform the suitability of forest harvesting as a source-water-protection strategy.
This study presents the results from an intensive field sampling program that characterized the water chemistry of multiple end members contributing to stream flow in legacy clear-cut and forested control catchments on the Canadian Shield. Increased knowledge of these runoff-generation processes is critical to explaining water quality responses to landscape disturbance. The specific objectives of this study were to (1) assess the primary end member contributions to stream flow in a legacy clear-cut (>20 years post harvest) and unharvested (control) catchments using end member mixing analysis and (2) assess temporal variability in end member contributions to stream flow during spring snowmelt, summer storm events, and fall wet-up.

Study Site
The Turkey Lakes Watershed (TLW; 47 • 03 N; 84 • 25 W) is an experimental watershed located in the Great Lakes in the St. Lawrence Forest Region of central Ontario (Figure 1). The watershed drains an area of 10.5 km 2 , and the relief is~300 m [41,43]. From 1980 to 2017, the mean annual temperature and precipitation was 4.5 • C and 1203 mm, respectively, with~35% of annual precipitation falling as snow [41,43]. The bedrock consists of Precambrian silicate greenstone (i.e., metamorphic basalt), with outcroppings of felsic igneous rock [41,56]. Overlying the bed rock is a two-component glacial till consisting of a compacted sandy basal till underlying a silt-loam ablation till [42,43]. Soils are classified as ortho humic-ferric podzols (spodosols), with well-defined L and F horizons typically 0.05 m thick [57]. The forest comprises predominantly uneven-aged, shade-tolerant, mature to over-mature hardwoods [56]. Sugar maple (Acer saccharum Marsh) comprises 90% of tree species, with yellow birch (Betula alleghaniensis Britton) and assorted conifers representing 9% and 1% of tree species, respectively [43,56].

Experimental Design
This study was conducted in two headwater sub-catchments within TLW. Catchment 31 (C31) has an area of 4.9 ha, total relief of 80 m, and a wetland area of ~2.0%. Catchment 32 (C32) has an area of 6.6 ha, total relief of 107 m, and a wetland area of ~1.5% [13]. In the fall of 1997, a silvicultural harvesting experiment was conducted in C31 to evaluate the effects of forest harvesting on soil productivity, stand recovery, biodiversity, and water quantity and quality [10]. C31 was treated with a diameter-limited or clear-cut harvesting strategy, which involved the felling and removal of all trees with a diameter at breast height (DBH) ≥ 20 cm. All trees with a DBH between 10 and 20 cm were felled and left on site, while trees with a DBH ≤ 10 cm were left standing [10,23]. Harvesting resulted in a 78% and 76% reduction in the stocking and basal area, respectively [36]. It should be noted that clear-cut harvesting is not a recommended harvesting strategy for Ontario's shadetolerant hardwood forests [58] but was selected to represent an extreme disturbance event [59]. C32 was not harvested to act as control because of its proximity and similar biophysical characteristics to C31. One notable difference between the two catchments is the position of the small, cryptic-treed wetlands (i.e., forested swamps) relative to the catchment outlet [60]. In C31, the wetland is in the upper parts of the catchment near the stream initiation point, whereas in C32, the wetland is located directly above the catchment outlet ( Figure 1).

Flow and Stream Water Chemistry
C31 and C32 were instrumented with 90° v-notch weirs, stilling wells, and Steven's Smart PT SDI-12™ pressure and temperature transducers to estimate stream flow at the basin outlet every 10 min. Stage-discharge rating curves were used to calculate instantaneous discharge. In addition, level loggers were placed in the upper parts of the stream

Experimental Design
This study was conducted in two headwater sub-catchments within TLW. Catchment 31 (C31) has an area of 4.9 ha, total relief of 80 m, and a wetland area of~2.0%. Catchment 32 (C32) has an area of 6.6 ha, total relief of 107 m, and a wetland area of~1.5% [13]. In the fall of 1997, a silvicultural harvesting experiment was conducted in C31 to evaluate the effects of forest harvesting on soil productivity, stand recovery, biodiversity, and water quantity and quality [10]. C31 was treated with a diameter-limited or clear-cut harvesting strategy, which involved the felling and removal of all trees with a diameter at breast height (DBH) ≥ 20 cm. All trees with a DBH between 10 and 20 cm were felled and left on site, while trees with a DBH ≤ 10 cm were left standing [10,23]. Harvesting resulted in a 78% and 76% reduction in the stocking and basal area, respectively [36]. It should be noted that clear-cut harvesting is not a recommended harvesting strategy for Ontario's shade-tolerant hardwood forests [58] but was selected to represent an extreme disturbance event [59]. C32 was not harvested to act as control because of its proximity and similar biophysical characteristics to C31. One notable difference between the two catchments is the position of the small, cryptic-treed wetlands (i.e., forested swamps) relative to the catchment outlet [60]. In C31, the wetland is in the upper parts of the catchment near the stream initiation point, whereas in C32, the wetland is located directly above the catchment outlet ( Figure 1).

Flow and Stream Water Chemistry
C31 and C32 were instrumented with 90 • v-notch weirs, stilling wells, and Steven's Smart PT SDI-12™ pressure and temperature transducers to estimate stream flow at the basin outlet every 10 min. Stage-discharge rating curves were used to calculate instanta-neous discharge. In addition, level loggers were placed in the upper parts of the stream channel near stream-initiation points. In C31, the stream originates from a small wetland referred to as C31 wetland, and in C32, the stream starts at the bottom of a hillslope referred to as the C32 initiation point. At these sampling locations, only the stream stage was recorded, as creating a stage-discharge rating curve was not possible due to the stream channel size and logistical constraints. Stream water samples were collected at both outlets and the two stream-initiation points using ISCO 6700 automated water samplers and stored in acid-washed, triple-rinsed 1 L polypropylene sampling bottles. ISCO samplers were programmed to collect daily composite (4 times a day) water samples during snowmelt and fall sampling periods. Additionally, the rising and falling limbs of select summer storms were sampled with ISCO samplers programmed to collect a sample every 0.5 to 2 h depending on the predicted size and length of the storm. Five storms were sampled during the summer. Finally, monthly water samples were collected by hand with 500 mL HDPE sampling bottles at all sampling locations to characterize baseflow conditions. All sampling bottles were acid-washed and triple-rinsed before being used in the field.

Groundwater
The ablation and basal till groundwater end member chemistry was characterized by collecting water samples from piezometer nests [37,38]. Piezometer nests consisted of two drive-point piezometers (Solinst Canada Ltd., Georgetown, ON, Canada) with 10 or 25 cm screens and an inside diameter (i.d.) of 2 or 4 cm, which were driven into the ablation or basal till [38]. Each catchment was instrumented with a total of 12 piezometer nests in two transects located in the upper and lower parts of the catchments. Each transect ran perpendicular to the stream channel, with three piezometer nests located on each side of the channel. Average depths of the piezometers were 0.45 ± 0.05 m and 0.40 ± 0.03 m in the ablation till and 1.09 ± 0.19 m and 0.93 ± 0.09 m in the basal till in C31 and C32, respectively [38]. In addition, wetland groundwater was characterized by a single transect of 4 cm i.d. piezometers installed to an average depth of 0.39 ± 0.03 m and 0.37 ± 0.03 m into the wetlands of C31 and C32, respectively. Groundwater samples collected using a peristaltic pump were obtained weekly during the snowmelt and fall sampling periods and during each storm event. Groundwater samples were also collected once per month during the summer to characterize groundwater chemistry during summer baseflow. Piezometers were purged for five minutes or until dry and then sampled the next day to allow the basal till piezometers sufficient time to recharge [37].

Soil Water
Soil pits with two zero-tension lysimeters were installed to collect soil water flowing through the soil profile next to each piezometer nest (excluding wetlands). Soil pits were dug by hand to an average depth of 0.65 ± 0.10 m and 0.64 ± 0.10 m in C31 and C32, respectively. Lysimeters were installed directly below the LFH layer and 20 cm below the LFH mineral soil interface. Lysimeter construction consisted of PVC pipe (15.23 cm i.d.), PVC sheeting, and a plastic screen mesh over the top to create a contact point between the lysimeter and soil/LFH layer. Additionally, the mineral soil lysimeters had glass wool, and a mixture of soil and deionized water was placed on top of the plastic screen mesh to strengthen the contact between the lysimeter and the soil matrix. Water drained from the lysimeters into acid-washed, triple-rinsed 1 L HDPE plastic sampling bottles. Sample bottles were placed in a plastic tote within each soil pit to protect them from the elements. All lysimeters were installed in the fall of 2020 to allow six months to equilibrate with the soil profile before sample collection in the spring of 2021. Two of the soil pits were completely saturated, preventing the installation of the lysimeters. As these soil pits were in small depressions, a drive-point piezometer was installed to a depth of 20 cm approximately one meter upslope of the soil pit to sample water flowing through the upper parts of the soil profile. Soil water samples were collected once a week during the snowmelt and fall sampling period and after each storm event. Fresh sampling bottles were placed in the soil pits before each storm to ensure any sample collected was representative of each storm. Lysimeter samples were also collected once a month in the summer corresponding with the groundwater sampling described above.

Throughfall
A throughfall collector was installed in each catchment to characterize precipitation chemistry. These collectors were constructed from plastic funnels with a mesh zip tied over the top to prevent clogging from leaf litter. The funnel was placed on a steel post next to a soil pit with a sample line running from the funnel to a 1 L HDPE sampling bottle within the tote containing the lysimeter sample bottles. Throughfall water samples were collected under the same sampling schedule as the lysimeters described above.

Laboratory Methods
All water samples were analyzed at the accredited Water Chemistry Laboratory in the Great Lakes Forestry Centre of the Canadian Forest Service in Sault Ste. Marie, Ontario, Canada. Water samples were analyzed for SO 4 2− , Ca 2+ , K + , Mg 2+ , Na − , Al 3+ , Ba 2+ , Sr 2+ , and DOC using standardized procedures and quality control methods [61]. Water chemistry tracers were selected based on tracers used in previous EMMA studies [48,49,53,54]. Due to the COVID-19 pandemic, there were several closures of the water chemistry laboratory, requiring sample freezing and storage for up to four months before processing and analysis could be completed.

End Member Mixing Analysis Model
End members used in EMMA model were selected based on physical characteristics of the catchments, existing monitoring equipment, and a prior understanding of the hydrological processes within TLW [23,37,38,43,60]. The specific end members selected included basal till groundwater, ablation till groundwater, wetland groundwater, mineral soil percolate, LFH percolate, and throughfall. Stream and end member sampling occurred during four different flow conditions defined as snowmelt, summer, summer storms, and fall.

Tracer Selection
Assumptions required for EMMA include the following: (i) tracer concentrations within end members need to be constant in time; (ii) tracers must be sufficiently different between end members; and (iii) tracers need to exhibit conservative behavior [45,46,48,49]. For the first assumption, end member tracers were grouped by season to remove any seasonal effects. The second assumption was tested using the Kruskal-Wallis H-test and tracer variability ratio (TVR). The Kruskal-Wallis H-test confirms for each tracer that at least one end member is significantly different (p ≤ 0.05) from another [48]. The TVR ratio, typically used in sediment source-fingerprinting studies, was used to determine if the difference in tracer concentrations between two end members was larger than the variation within either end member [62]. TVR was calculated using Equation (1) for each tracer and each end member group pairing: whereX max andX min are the maximum and minimum median tracer concentrations of the two potential end members, and CV is the coefficient of variation of either end member [62][63][64]. A TVR ratio ≥ 2 was required for the tracer to be included in the EMMA, as recommended by [62,63]. The third assumption was evaluated using bivariate plots that are commonly used to assess conservative mixing in EMMA studies [47,48,50,63]. Linear relationships between tracers suggest conservative behavior, as nonlinearity may be evidence of chemical reactions occuring between solutes [46]. A Pearson's r correlation coefficient ≥ 0.5 (p < 0.01) was used as the cutoff to represent a linear relationship between two tracers [48], and a tracer that had a linear relationship with at least one other tracer was deemed acceptable. All tracers that passed the critical assumptions unperdinning EMMA were considered appropriate.

End Member Mixing Analysis Model
The EMMA model was created following the methods of Christopherson and Hooper [46] and Hooper [50]. For detailed, step-by-step instructions, the readers can refer to any of the following [48,51,53,65]. The mixing sub-space (or U space) was created for each sampling location under the different flow conditions to create a two-dimensional (2D) mixing sub-space [46]. Two dimensions were selected for ease of visualization [63]. Stream water chemistry (n samples × p tracers) was standardized and centered (subtracted by the mean and divided by the standard deviation) to create equal weighting between tracers. PCA was conducted on the correlation matrix of the standardized and centered data using the prcomp function in R [66]. PCA analysis was completed on all tracers, and every combination of four up to p tracers was computed, and the PCA model that explained the most variation in two dimensions was retained [53,54]. Median ± median absolute deviation (MAD) end member tracer concentrations were then standardized and centered around the stream water chemistry and projected into the mixing subspace defined by the stream water chemistry. If all end members were characterized correctly, and conservative mixing did occur, then end members were bound within the stream water chemistry.
The diagnostic tool described by Hooper [49] was used to confirm that a twodimensional mixing subspace was appropriate. Residual analysis between the original stream water chemistry and the orthogonal projections for each solute was completed, with a random pattern between the residuals and the original stream water chemistry describing a good mixing subspace. If a random pattern was not observed, then additional eigenvectors (dimensions) were retained until a random pattern was observed [48,50,67]. Additionally, the fit between the orthogonal projections and observed stream water chemistry was evaluated using the relative root-mean-square error (RRMSE) for each solute and dimension: where n represents the number of stream water samples, X p represents the median concentration of solute p,X ip is the orthogonal projection from the mixing space, and X ip is the observed concentration of solute p in sample i. Tracers with a RRMSE < 15% were deemed appropriate for EMMA [47].

Hydrometeorological Characterization
Stream flow sampling during snowmelt occurred from 25 March 2021 to 9 May 2021 and captured the latter half of snowmelt, which typically starts at the beginning of March in these study catchments [41]. Total precipitation during this period was 117.5 mm, and runoff in the clear-cut and control catchments was 140.1 mm and 151.6 mm, respectively. Maximum daily runoff for the clear-cut and control catchments was 16.6 mm day −1 and 16.6 mm day −1 , respectively (Table 1). Fall sampling occurred from 9 October 2021 to 30 October 2021, with a total precipitation of 97.9 mm and total stream flow of 18.5 mm and 27.4 mm in the clear-cut and control catchments, respectively. Peak daily stream flow was 5.1 mm day −1 in the clear-cut and 2.7 mm day −1 in the control catchments ( Table 1). The number of zero-flow days observed during the field sampling period at the four sampling locations was highest at the C32 initiation point, with 124 zero-flow days ( Table 2). The clear-cut and clear-cut wetland experienced 48 and 65 zero-flow days, respectively, while the control catchment had no zero-flow days at its outflow (Table 2).   (Table 3). Total precipitation for storm 4 was 14.1 mm, while total stream flow was 2.8 mm for the clear-cut and 1.9 mm for the control. Peak instantaneous stream flow was 0.7 L s −1 ha −1 and 0.2 L s −1 ha −1 for the clear-cut and control catchments, respectively (Table 3). Antecedent 2-, 7-, and 14-day precipitation metrics were 25.7 mm, 26.0 mm, and 41.8 mm for storm 1 and 14.1 mm, 22.2 mm, and 43.9 mm for storm 4 (Table 3).

Tracer Selection
The number of samples needed for each end member during the four flow conditions was evaluated to ensure sufficient data were available for end member characterization.
Few observations were available for the ablation and basal till end members during the summer storm sampling periods (Table 4) because the groundwater table was often below the ablation and basal till piezometers. Accordingly, the summer and summer storm data were combined to create the end members for this period. It was hypothesized that there would be no difference in tracer concentration between these two sampling periods, as they both occurred during the dry summer months. It should be noted that even after combining the data, only one datapoint was available for the ablation and basal tills in the control catchment. This is not sufficient to characterize these end members, but they were still included in EMMA, as the variation around the groundwater end members is small (Figure 2). To be consistent with the changes made for the ablation and basal till end members, all end member data for the summer and summer storms sampling periods were combined. Significant differences (Kruskal-Wallis H-test p ≤ 0.05) were found between all end members for all tracers in both the clear-cut and control catchments (Figure 2). TVR ratios were >2 for all end member comparisons for all tracers for both clear-cut and control catchments. A bivariate plot of stream water chemistry in the clear-cut during snowmelt is shown in Figure 3 and highlights that each tracer behaves conservatively, as each has a Pearson's r correlation coefficient ≥ 0.5 [48]. All bivariate plots can be seen in Figures S1-S11. A summary of the tracer selection processes is shown in Table 5.

End Member and Stream Water Characterization
Solute concentrations ordered by throughfall, LFH percolate, wetland groundwater, mineral soil percolate, ablation till, basal till, and stream water ( Figure 2) indicate that concentrations either increased or decreased sequentially from throughfall to the basal till groundwater end members. Generally, concentrations of Ca 2+ , Mg 2+ , Na − , SO4 2− , and Sr 2+ were highest in the ablation and basal till groundwater end members and lowest in the throughfall and LFH percolate end members ( Figure 2). Conversely, Al 3+ , K + , and DOC were elevated in the throughfall and LFH end members and lower in the ablation and basal till groundwater end members ( Figure 2). One exception to the above patterns was Ba 2+ , which showed the highest concentrations in the wetland groundwater and mineral soil percolate end members and the lowest concentrations in the throughfall end member (Figure 2). There were few differences in these patterns between the clear-cut and control catchments and the snowmelt, summer, and fall flow conditions. Tracer concentrations that were highest in the groundwater end members (Ca 2+ , Mg 2+ , Na − , SO4 2− , and Sr 2+ ) increased throughout snowmelt in both the control and clear-cut catchments ( Figure 4) while showing episodic decreases with increasing flow rates. During the fall wet-up, stream concentrations of Ca 2+ , Mg 2+ , Na − , SO4 2− , and Sr 2+ decreased as flow increased and then showed no change as flow stabilized. During the summer storms, Ca 2+ , Mg 2+ , Na − , SO4 2− , and Sr 2+ were inversely related to increasing stream flow ( Figure 5). The opposite pattern was observed for Al 3+ , K + , and DOC, which were highest in the throughfall and LFH end members and increased with increasing flow during snowmelt, fall, and summer storms (Figures 4 and 5). No pattern was observed for Ba 2+ in stream water for any sampling period (Figures 4 and 5).

End Member and Stream Water Characterization
Solute concentrations ordered by throughfall, LFH percolate, wetland groundwater, mineral soil percolate, ablation till, basal till, and stream water ( Figure 2) indicate that concentrations either increased or decreased sequentially from throughfall to the basal till groundwater end members. Generally, concentrations of Ca 2+ , Mg 2+ , Na − , SO 4 2− , and Sr 2+ were highest in the ablation and basal till groundwater end members and lowest in the throughfall and LFH percolate end members ( Figure 2). Conversely, Al 3+ , K + , and DOC were elevated in the throughfall and LFH end members and lower in the ablation and basal till groundwater end members ( Figure 2). One exception to the above patterns was Ba 2+ , which showed the highest concentrations in the wetland groundwater and mineral soil percolate end members and the lowest concentrations in the throughfall end member (Figure 2). There were few differences in these patterns between the clear-cut and control catchments and the snowmelt, summer, and fall flow conditions. Tracer concentrations that were highest in the groundwater end members (Ca 2+ , Mg 2+ , Na − , SO 4 2− , and Sr 2+ ) increased throughout snowmelt in both the control and clear-cut catchments ( Figure 4) while showing episodic decreases with increasing flow rates. During the fall wet-up, stream concentrations of Ca 2+ , Mg 2+ , Na − , SO 4 2− , and Sr 2+ decreased as flow increased and then showed no change as flow stabilized. During the summer storms, Ca 2+ , Mg 2+ , Na − , SO 4 2− , and Sr 2+ were inversely related to increasing stream flow ( Figure 5). The opposite pattern was observed for Al 3+ , K + , and DOC, which were highest in the throughfall and LFH end members and increased with increasing flow during snowmelt, fall, and summer storms (Figures 4 and 5). No pattern was observed for Ba 2+ in stream water for any sampling period (Figures 4 and 5).

EMMA Results
Stream water chemistry was not bounded by end members in the dataset (mixing subspace), and therefore, hydrograph separation was not possible. However, prominent end members contributing to the stream water mixture were identified based on their proximity to stream water values as described by Spencer et al. [63].

Snowmelt
During snowmelt, two principal components (PCs) explained 97.4%, 98.2%, and 95.5% of the variance for the clear-cut, control, and clear-cut wetland, respectively ( Figure 6). The ablation and basal till groundwater end members were most like stream water chemistry for all three sampling locations, with the wetland end member being the next-closest for the clear-cut and clear-cut wetland. No wetland groundwater data were available in the control catchment during snowmelt. Throughfall and LFH percolate end members were most dissimilar from the stream water for all three end members. Higher flow rates in the clear-cut and control catchments trended towards the LFH percolate and mineral soil percolate end members. No obvious pattern was observed for increasing stream stage in the clear-cut wetland, and loadings trended away from the end members. Exceptions include K + in the clear-cut and control trending towards the mineral soil and LFH percolate end members and DOC in the clear-cut wetland trending towards wetland groundwater, mineral soil, and LFH percolate in PC2.

Summer Storms
In the clear-cut, control, and clear-cut wetland, two PCs explained 97.8%, 96.6%, and 97.3% of the variance during storm 1 (Figure 7). Stream water values were closest to the ablation and basal till end members in the clear-cut, control, and clear-cut wetland. Additionally, the wetland groundwater and mineral soil water percolate were close to the stream water values in the control and mineral soil percolate in the clear-cut wetland. Increasing stream flow trended towards the wetland groundwater, throughfall, mineral soil, and LFH percolate in the clear-cut and clear-cut wetland. In the clear-cut catchment, Al 3+ trended towards the wetland groundwater, mineral soil, and LFH percolate end member, while Mg 2+ and Ca 2+ trended towards the ablation and basal till end members. In the control catchment, Mg 2+ and Ca 2+ trended towards the wetland groundwater end member, while K + trended towards the throughfall, mineral soil, and LFH percolate end members. For the clear-cut wetland, SO 4 2− trended towards the ablation till end member, while Mg 2+ and Ca 2+ trended towards the basal till end member.
For storm 4, two PCs explained 98.3%, 97.3%, and 98.2% of the variance in the clearcut, control, and clear-cut wetland, respectively (Figure 7). Stream water chemistry was closest to the ablation and basal till end members in the clear-cut and clear-cut wetland, while wetland groundwater and ablation till end members were closest in the control. Increasing stream flow trended towards the mineral soil and LFH percolate in the clear-cut and towards all end members in the control. DOC trended towards the mineral soil and LFH percolate in the clear-cut, while all loadings trended away from end members in the control. DOC and Al 3+ trended towards the wetland groundwater and LFH percolate in the clear-cut wetland.    . EMMA models for the clear-cut, control, and clear-cut wetland sampling locations during storm 1 and storm 4 using the first two principal components that explain the most variation in stream water chemistry. Grey circles show stream water points, with circle size denoting either stream flow (L s −1 ha −1 ) for the clear-cut and control catchment or head (cm) for the clear-cut wetland. Colored points and error bars denote the median ± MAD end member tracer concentrations that were projected into the mixing subspace defined by the stream water chemistry. Arrows denote the loadings from the principal component analysis (PCA) for each tracer used in the EMMA model.

Fall
Two PCs explained 99.1%, 99.4%, and 99.2% of the variance in the clear-cut, control, and clear-cut wetland, respectively ( Figure 6). The ablation and basal till groundwater end members were closest to the stream water values at the clear-cut weir and clear-cut wetland, while the mineral soil percolate end member was closest to the stream water values in the control catchment. Higher flow rates trended towards the LFH percolate and wetland groundwater in the clear-cut and control catchments, respectively. Most loadings trended away from the end members. Exceptions include K + moving towards wetland groundwater in the control and SO 4 2− towards all end members in PC1 for the clear-cut wetland.

Residuals Analysis
Residuals analysis (RRMSE < 15%) confirmed that a two-dimensional mixing space was sufficient to explain the stream water chemistry for most sample locations for a range of flow conditions (Figure 8). Exceptions included the clear-cut wetland during storm 4, as DOC was above the 15% cut-off in the EMMA model ( Figure 8). If stream water chemistry was properly bound by the end members, then a three-dimensional mixing space (four end members) would be required to explain the stream water chemistry in the clear-cut wetland during storm 4. wetland, while the mineral soil percolate end member was closest to the stream water values in the control catchment. Higher flow rates trended towards the LFH percolate and wetland groundwater in the clear-cut and control catchments, respectively. Most loadings trended away from the end members. Exceptions include K + moving towards wetland groundwater in the control and SO4 2− towards all end members in PC1 for the clear-cut wetland.

Residuals Analysis
Residuals analysis (RRMSE < 15%) confirmed that a two-dimensional mixing space was sufficient to explain the stream water chemistry for most sample locations for a range of flow conditions (Figure 8). Exceptions included the clear-cut wetland during storm 4, as DOC was above the 15% cut-off in the EMMA model ( Figure 8). If stream water chemistry was properly bound by the end members, then a three-dimensional mixing space (four end members) would be required to explain the stream water chemistry in the clearcut wetland during storm 4. Figure 8. Relative root-mean-square error (RRMSE, %) for the tracers used in the EMMA models. Horizontal line represents the 15% cut-off threshold recommended by James and Roulet [46]. Color represents the RRMSE explained by each principal component or dimension of the mixing subspace.

Validity of EMMA at TLW
Three end members are required to explain stream water chemistry in the catchments (Figure 8) [50], and the data provided herein are comparable to other studies [48,52]. However, hydrograph separation in the present study was not completed because the stream water was not bound by end members within the mixing subspace. The only scenario where hydrograph separation could have been calculated was for the clear-cut wetland and the control catchment during storm 1 (Figure 7). However, this analysis still provides key insight into hydrologic processes that control runoff generation within these Horizontal line represents the 15% cut-off threshold recommended by James and Roulet [47]. Color represents the RRMSE explained by each principal component or dimension of the mixing subspace.

Validity of EMMA at TLW
Three end members are required to explain stream water chemistry in the catchments (Figure 8) [50], and the data provided herein are comparable to other studies [48,52]. However, hydrograph separation in the present study was not completed because the stream water was not bound by end members within the mixing subspace. The only scenario where hydrograph separation could have been calculated was for the clear-cut wetland and the control catchment during storm 1 (Figure 7). However, this analysis still provides key insight into hydrologic processes that control runoff generation within these catchments and forms the conceptual framework to explain changes in stream chemistry. Studies have reported that while end members were not bound by stream water chemistry, they were previously used to infer runoff-generating processes in the Rocky Mountains of Alberta [63] and Canadian Shield of Quebec [48]. The lack of bounding of stream water by end member tracer concentrations could be caused by a missing end member that contributes to stream water [48,63]. Figures 2 and 4 show that concentrations of Ca 2+ , Mg 2+ , Na − , SO 4 2− , and Sr 2+ within stream water lay outside the range observed in all the end members, which suggests that another end member is missing from this analysis. Missing end members could include a deeper groundwater source or a specific hillslope position that contribute to preferential inflows, which was not sampled in this study. Ali et al. [48] discussed the important role that aspect and specific hillslope positions, such as gullies, played in contributing to stream flow in the hardwood forests of the Hermine catchment at the Station de Biologie des Laurentides, Quebec, Canada, a study area having similar climate, geology, soils, and vegetation to the TLW. Additionally, specific hillslope positions or preferential inflows have been shown to influence stream water quality parameters such as DOC [68] and temperature [69]. Notably, the sampling design employed in this study may have precluded measurement of the complete range of water quality parameters for soil water and preferential surface flow pathways that contribute end members to stream flow.
It has been suggested that EMMA is not appropriate for hydrograph separation because spatial and temporal variation in solute end member chemistry with changing antecedent moisture conditions violates the assumptions of this procedure [70]. In the present study, end member chemistry was split by season to minimize these changes; however, the variability in solute chemistry (Figures 6 and 7) and shifting median end member solute concentrations ( Figure 4) suggest these assumptions were invalidated [48,70]. Accordingly, this provides further support that hydrograph separation using the mass balance equations should be avoided. Notably, this observation combined with the laborintensive sampling regime required for EMMA may lead to the perception that it is not a cost-effective method for discerning hydrologic processes. However, using the proximity of end members to stream water chemistry within the mixing subspace (Figures 6 and 7) and the hydro-chemical response of stream water (Figures 4 and 5), inferences on hillslope channel connectivity and runoff-generation processes were still obtained. Additionally, this study highlights the importance of utilizing appropriate and thus different EMMA models (i.e., different tracers or end members) to represent different hydrologic scenarios, as shown by Ali et al. [48].

Insights into Hydrologic Processes
Ablation and basal till groundwater chemistries were the end members most similar to stream water in the mixing subspace, highlighting their role as primary sources contributing to stream water. This pattern was consistent across all three sampling locations with few seasonal differences (Figures 6 and 7). This agrees with previous work at TLW, which showed that during snowmelt, a perched water table over the basal till routed water through the ablation till before reaching the stream channel [42,43]. Stream water chemistry at high flows trended towards the wetland groundwater, mineral soil water, and LFH percolate (Figures 6 and 7), suggesting that the upper portions of the soil profile are only hydrologically connected during higher-flow conditions. Figure 4 supports this, with Al 3+ , K + , and DOC concentrations, which are highest in the LFH percolate, mineral soil water, and wetland groundwater end members, peaking in stream water at higher flows. Similar studies have shown that groundwater is often the dominant contributing end member to stream water in forested headwater catchments [52,53,65], and only under higher-flow conditions do the upper parts of the soil profile and hillslope become hydrologically connected [48,71].
Few seasonal differences in end member contributions to stream flow were observed and were rarely consistent among the three sampling locations. Exceptions were observed with the wetland groundwater and mineral soil percolate end members being closest to stream water values for the control catchment during the fall and summer storms (Figures 6 and 7). The lower contributions from the ablation and basal till may be due to the lack of data from these end members in the control catchment during summer. Given that few groundwater samples were collected during storm events, and a stream response to storms was still observed, it could be hypothesized that these end members were not major contributors to stream water, as there was rarely enough ground water to collect a sample. The location of the piezometers may have also missed critical subsurface flow paths such as those described by Ali et al. [48]. Additionally, the proximity of the wetland to the catchment outlet explains the prominent contributions of wetland groundwater to stream flow in the control catchment.

Wetland Contributions to Stream Flow
While groundwater contributions were the primary source of water within the clear-cut and control catchments, wetland groundwater contributions were an important secondary source, especially during summer storms (Figures 6 and 7). This is especially apparent in the control catchment, where the wetland is directly upstream of the catchment outlet. During the summer storms, wetland groundwater solute chemistry was closest to stream water in the mixing subspace in the control catchment, while the groundwater end members were closest in the clear-cut (Figure 7). The wetland signal within the clear-cut may be diluted by hillslope contributions, as the wetland is near the stream-initiation point, unlike its location directly above the catchment outlet in the control. Further support for this hypothesis is the number of zero-flow days observed in these catchments. Flow was sustained in the control catchment during the summer (no zero-flow days), while the clear-cut had 48 zero-flow days ( Table 2). While these differences may be caused by the impacts of legacy harvesting, they are more likely the result of differences in wetland position between the two study catchments. Studies have shown that small wetlands can act as sponges that fill up during wet periods and slowly drain, sustaining flow during dry periods [24,72]. Work exploring the effect of wetland cover on mean water-travel times supports this theory, as the relationship between travel times and wetland cover is dependent on the antecedent moisture conditions [23,73]. During wet time periods, wetlands rapidly transport water laterally through the upper parts of the soil profile due to low water-storage capacity and high hydraulic conductivities [23]. During dry periods, wetlands may act as a small reservoir collecting water from surrounding hillslopes and slowly draining into the stream channel [23,73]. Notably, in the present study, stream flow within the control catchment was mostly sourced by wetland groundwater. The proximity of the wetland to the control catchment outlet maintained flow throughout the summer.

Legacy Harvesting Effects on End Member Contributions to Stream Flow
Despite the differences in wetland groundwater contributions to stream flow between the control and clear-cut catchments, there are few differences between the clear-cut and control catchments, suggesting legacy forest harvesting has limited impact on the prominent end member contributions to stream flow. Previous work at TLW found that forest harvesting increased the proportion of water moving through the upper parts of the soil profile based on K:SiO 2 molar ratios four years after harvesting [38]. If harvesting was still impacting these flow paths, it would be expected that the clear-cut mineral soil and LFH percolate end member contributions to stream water would be larger than in the control. As no differences were observed in end member contributions to stream flow between the control and clear-cut catchments, it can be assumed that the harvesting-induced flow path changes have recovered. As signs of flow path recovery have been observed 15 years after harvesting [10], it is reasonable to assume a complete recovery may have occurred by the time of this study. This suggests that forest harvesting may be a suitable source-waterprotection strategy for this landscape, as there are limited long-term harvesting impacts on runoff-generation processes that most affect stream water quality. However, recent work conducted at TLW has shown how harvesting impacts on water yield in an adjacent catchment that was subject to a selection harvest treatment may have persisted until 2020 [10,74]. This in combination with the high variability observed in end member concentrations suggests the results from this study should be interpreted cautiously. However, the lack of change in end member contributions to stream flow does not necessarily mean there would be no change in stream water yield. Future work should continue to explore the impacts of forest harvesting on these runoff-generation processes to further inform the impacts on both water quantity and quality.

Implications for Water Quality and Treatability
Previous work has provided insights to the behavior of numerous water quality parameters and their response to harvesting at TLW [13,61]. Solutes such as K + and DOC were still elevated in the study catchments 21 years after harvesting; these responses were attributed to long-term changes in solute demand/availability and hillslope channel connectivity [13]. As both solutes are generally associated with shallow flow paths, it follows that mineral soil and LFH end members may contribute more to stream flow in harvested catchments. However, the results from EMMA did not support this hypothesis and may suggest that changes in solute availability on the landscape are driving the observed harvesting response. This shows the importance of quantifying not only hillslope channel connectivity but also nutrient and solute availability when evaluating the suitability of landscape management decisions such as forest harvesting as source-water-protection strategies.
Previous studies have identified wetlands as a critical source of numerous water quality parameters at TLW [23,60,75]. Solutes such as DOC, dissolved organic nitrogen (DON), total phosphorus (TP), and total dissolved phosphorus (TDP) were all elevated in catchments with large wetland fractions at TLW [60,75]. The unique biogeochemical processes occurring in these wetlands heavily influence the solubility, lability, and mobility of these nutrients [60,75]. The results from this study clearly show the important influence wetlands play in runoff generation within these headwater catchments, providing further evidence of their influence on both water quantity and quality. Wetlands appear to have a disproportionately larger influence on stream water quantity and quality relative to other landscape features such as hillslopes. Therefore, special care should be taken by land managers to ensure wetlands protection given their sensitivity to landscape disturbance [24].

Conclusions
End member mixing analysis was used to evaluate the dominant sources of hillslope runoff that contribute to the generation of stream flow in a legacy clear-cut and control catchment at TLW. Tracer concentrations either increased (Ca 2+ , Mg 2+ , Na − , SO 4 2− , and Sr 2+ ) or decreased (Al 3+ , K + , and DOC) sequentially from throughfall to the basal till groundwater end members. Stream water was not bounded by end members in the mixing subspace, making hydrograph separation impossible using this dataset. Despite this limitation, results suggest that ablation and basal till groundwater end members were the major contributors to stream flow in both the legacy clear-cut and control catchments during snowmelt and fall. Groundwater originating from wetlands within the stream network was also identified as a prominent source of stream flow in both catchments. During the summer storms, wetland groundwater was the dominant source of stream flow to the control catchment. However, wetland position in the landscape and the degree of hydrologic connectivity to the stream channel influenced the contributions of wetland groundwater to stream flow. There were few differences between end member contributions to stream flow in the legacy clear-cut and control catchments suggesting hydrologic recovery had largely occurred 24 years post harvesting. Future studies should focus on identifying critical hillslope positions or preferential inflows that likely represent large contributions to stream flow within these catchments. Identifying these critical positions may increase the chances of hydrograph separation using EMMA and provide more quantitative insight into the impacts of changing antecedent moisture conditions and landscape disturbances on runoff-generation processes. Overall, the results from this study suggest that legacy forest harvesting has a relatively small impact on hydrological processes. Accordingly, forest harvesting may be a suitable source-water-protection strategy for Canadian Shield catchments. Additionally, this study demonstrates the utility of EMMA to evaluate key drivers that influence the relationship between hydrologic processes and water quality. Such information is critical for land managers and drinking water providers to develop and evaluate source-water-protection strategies and drinking-water safety plans and inform landscape management decisions.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/w15152825/s1, Figure S1: Bivariate plot for the clear-cut during the fall; Figure S2: Bivariate plot for the clear-cut during storm 1; Figure S3: Bivariate plot for the clear-cut during storm 4; Figure S4: Bivariate plot for the control during snowmelt; Figure S5: Bivariate plot for the control during the fall; Figure S6: Bivariate plot for the control during storm 1; Figure S7: Bivariate plot for the control during storm 4; Figure S8: Bivariate plot for the clear-cut wetland during snowmelt; Figure