Drivers of Dissolved Organic Carbon Mobilization From Forested Headwater Catchments: A Multi Scaled Approach

Understanding and predicting catchment responses to a regional disturbance is difficult because catchments are spatially heterogeneous systems that exhibit unique moderating characteristics. Changes in precipitation composition in the Northeastern U.S. is one prominent example, where reduction in wet and dry deposition is hypothesized to have caused increased dissolved organic carbon (DOC) export from many northern hemisphere forested catchments; however, findings from different locations contradict each other. Using shifts in acid deposition as a test case, we illustrate an iterative “process and pattern” approach to investigate the role of catchment characteristics in modulating the steam DOC response. We use a novel dataset that integrates regional and catchment-scale atmospheric deposition data, catchment characteristics and co-located stream Q and stream chemistry data. We use these data to investigate opportunities and limitations of a pattern-to-process approach where we explore regional patterns of reduced acid deposition, catchment characteristics and stream DOC response and specific soil processes at select locations. For pattern investigation, we quantify long-term trends of flow-adjusted DOC concentrations in stream water, along with wet deposition trends in sulfate, for USGS headwater catchments using Seasonal Kendall tests and then compare trend results to catchment attributes. Our investigation of climatic, topographic, and hydrologic catchment attributes vs. directionality of DOC trends suggests soil depth and catchment connectivity as possible modulating factors for DOC concentrations. This informed our process-to-pattern investigation, in which we experimentally simulated increased and decreased acid deposition on soil cores from catchments of contrasting long-term DOC response [Sleepers River Research Watershed (SRRW) for long-term increases in DOC and the Susquehanna Shale Hills Critical Zone Observatory (SSHCZO) for long-term decreases in DOC]. SRRW soils generally released more DOC than SSHCZO soils and losses into recovery solutions were higher. Scanning electron microscope imaging indicates a significant DOC contribution from destabilizing soil aggregates mostly from hydrologically disconnected landscape positions. Results from this work illustrate the value of an iterative process and pattern approach to understand catchment-scale response to regional disturbance and suggest opportunities for further investigations.


INTRODUCTION
Disturbances such as land-cover transformation, amplification of biogeochemical flows, and climate disruption are triggering transitions in the Earth system that are unprecedented on human timescales (Steffen et al., 2018;Abbott et al., 2019). Potential and ongoing ecosystem state changes may impact billions of individuals (Abatzoglou and Williams, 2016;Van Loon et al., 2016;Dupas et al., 2019), highlighting the need to better understand the factors that determine how catchments respond to multiple disturbance types and in different contexts.
However, attributing catchment response to even a specific driver is challenging because the scale of disturbance often transcends the bounds of individual sites that, themselves, have heterogeneous spatial characteristics. While spatial and temporal patterns at larger scales can be extracted from relevant data (e.g., atmospheric deposition data), finding patterns that hold true across individual catchments remains a challenge that was already described decades ago as "uniqueness of place" for hydrological modeling (Beven, 2000). Research shows that a direct scaling of site-specific observations to a larger-scale trend (and vice versa) might not be possible and that scaling mismatches might be the norm rather than the exception (Levin, 1992;NSF, 2018). For example, Preston et al. (2011) used lab mesocosm, field monitoring and long-term stream data to investigate the impact of climate drivers on dissolved organic carbon (DOC) liberation, and found that variations in temperature and moisture impacted DOC concentrations differently for each scale of observation. While direct upscaling and downscaling might lead to loss of information and insight, an approach that does not seek direct links across scales but that investigates larger-scale pattern and site-specific process in an iterative fashion might provide some remedy (Sivapalan, 2005). In such an approach, regional patterns might point to a possible process (i.e., pattern-to-process), while process investigation might provide more insights with each iteration (i.e., process-to-pattern).
Acidification (and ensuing reduction in acid deposition) is a good candidate for iterative pattern and process investigation in the context of disturbances. Large variation in forested catchment signals including increases in stream DOC exports are attributed reductions in acid deposition, but the pattern is not consistent (Driscoll et al., 2003;Burns et al., 2006;De Wit et al., 2007;Monteith et al., 2007;Hruška et al., 2009;Sanclements et al., 2018). For example, some catchments have shown increases in DOC concentration despite continued acid deposition (Oni et al., 2013), while others showed reduced deposition without attendant increases in DOC (Löfgren and Zetterberg, 2011). Concurrent with the reduction in acid deposition, changes to the climate system have led to increasing precipitation, which is amplified by increasing frequency of extreme hydrological events such as heavy precipitation (Jentsch et al., 2007;Campbell et al., 2009;Arnone et al., 2011;Seneviratne et al., 2012). For example, in the case of catchments in the Northeastern US, increased heavy precipitation and reduced acid deposition could lead to increased exports of DOC, as shown in a long-term, paired catchment study (Sanclements et al., 2018) where one catchment received continuous acid treatment while the other was allowed to recover from acidification. The latter catchment exported significantly more DOC, indicating a strong effect of reduced acid deposition; however, the temporal variability in stream DOC concentrations for both catchments was high due to shifts in precipitation patterns. While these overlapping disturbances alone can lead to DOC patterns that are difficult to interpret, variations in catchment-specific characteristics likely drive differences in DOC export patterns in areas that receive similar acid loadings (Clark et al., 2010). The authors of the latter study note that seasonal, short-term inter-annual and longerterm temporal variation interact in complex ways, masking otherwise compatible patterns, and suggest careful process investigation. Lastly, variation in length of record and reporting methods (Eimers et al., 2008a) present a continuous challenge.
At the catchment scale, processes that lead to DOC liberation into streams are coupled and highly complex. Especially in forested headwater catchments, DOC is primarily controlled by soil organic matter (SOM) that supplies DOC when soils are flushed in response to precipitation and, in snow-dominated systems, during snowmelt (Boyer et al., 1996;Raymond and Saiers, 2010;Wilson et al., 2013;Moatar et al., 2017;Perdrial et al., 2018;Zarnetske et al., 2018;Zhi et al., 2020). Thus, DOC is often controlled by discharge (Q) and studies now increasingly use flow-adjusted DOC concentrations (DOC FA ) or fluxes in their assessment to test for DOC drivers other than Q (Burns et al., 2006). This hydrological control is also embedded in the broader context of spatial and temporal variability of catchment processes. Recent observations and reactive transport modeling support the shallow and deep hypothesis that high concentration DOC in stream under wet, high discharge conditions comes from shallow soils enriched with SOM, whereas low stream DOC under dry, low discharge conditions originates from deeper groundwater that is often depleted in DOC (Zhi et al., 2019). The concept of ecosystem control points is very useful in exemplifying this complexity and heterogeneity because variable importance of DOC supply vs. transport is acknowledged (Bernhardt et al., 2017). For example, catchment locations such as planar hillslopes accumulate materials until they are hydrologically connected and episodically flushed (export control points) and lead to high temporal variability in DOC exports in streams during events. In contrast, locations such as riparian zones and swales exhibit high degrees of biogeochemical activity and near continuous connectivity to the stream (permanent control points) and the ensuing impact on stream DOC is more continuous. Because microbial processing is regulated by soil temperature and moisture (Wen et al., 2020), landscape positions with favorable soil moisture conditions promote microbial DOC processing and losses to the atmosphere as CO 2 , thus reducing DOC concentrations in soil solution (Moyano et al., 2013). In turn, very wet conditions suppress aerobic DOC respiration leading to DOC accumulation available for flushing (Huang and Hall, 2017). Because SOM is an important supplier for DOC, SOM stability also needs to be considered in this context. SOM can be stabilized chemically, through adsorption on clay and silt particles (Six et al., 2002;Perdrial et al., 2010) or physically, through the occlusion within soil micro-aggregates (<250 µm) that limit microbial processing (Lin et al., 2006;Mikutta et al., 2006;Kleber and Johnson, 2010;Schmidt et al., 2011;Li et al., 2017). Recent research has indeed begun to link reduced acid deposition and subsequent changes in soil chemistry (ionic strength and pH) to changes in SOM stabilization in aggregates driving DOC loss at the catchment scale (Armfield et al., 2019;Cincotta et al., 2019). Because these conditions vary in scale, both in space (region, catchment, landscape position within a catchment) and time (decadal vs. seasonal), this complexity is difficult to detangle and a good test case for iterative pattern and process investigations.
The overarching goal of this study is to illustrate opportunities and challenges with the process and pattern approach using increased and decreased acid deposition as test case. We apply one iteration as starting point with the intent to critically evaluate the opportunities and limitations of each approach, as well as the integration of observations across these scales. For pattern investigation at the regional scale, we investigate linkages between reduced acid deposition and stream DOC response and the potential role of catchment characteristics as covariates in that response (objective 1). For this we apply regional atmospheric deposition data from the National Atmospheric Deposition Program (NADP), catchment biogeophysical and hydrometeorological attributes (Addor et al., 2017), and a newly curated dataset that includes co-located stream Q and stream chemistry data (Sterle et al., 2019). We use flow-adjusted stream chemistry data to remove the typically dominant Q control from temporal trends and to better isolate potential signals in response to reduced acid deposition in stream DOC concentrations (Helsel et al., 2020). We compute trend analyses for long-term data sets and compare stream water DOC trend directionality against sulfate (SO 2− 4 ) deposition trends as one proxy for shifts in atmospheric deposition, and we compare these trends to catchment attributes (Addor et al., 2017) to determine if common catchment scale characteristics might be associated with these trends. For process investigations, we test the effect of leaching solution composition on soil aggregate stability and DOC release experimentally (objective 2). For this, we use soils cores from several landscape positions and locations in catchments with opposite long-term DOC trends [Susquehanna Shale Hills Critical Zone Observatory (SSHCZO) and the Sleepers River Research Watershed (SRRW)] and critically evaluate opportunities and limitations of integrating these observations with large-scale patterns (Figure 1).

METHODS AND MATERIALS
Long-Term Trend Analyses at the Regional Scale

Site Selection and Data Collection
To characterize trends in DOC and Q, we utilized the recently-developed CAMELS-Chem dataset (Sterle et al., 2019). This dataset complements the existing CAMELS dataset with additional atmospheric deposition and water chemistry data for 671 catchments, 499 of which are >50% forest cover. The temporal record of stream chemistry varies depending on the individual catchment, with all records ranging between 1980 and 2019. In this study, we constrained the dataset to forested catchments, that is catchments with >50% forest cover, in the Northeastern United States. This was done to minimize land cover drivers from our analyses (Figure 2). Within the Northeast, we further constrained catchments according to their length of available stream DOC and Q records (at least five sequential years) and the density of the records (i.e., no gaps >1/3rd of the total record length) (Meals et al., 2011). Once we identified suitable catchments, we compiled all available DOC concentration and Q records from said catchments for trend and statistical analyses.
To characterize trends in wet deposition of SO 2− 4 , we utilized the National Atmospheric Deposition Program (NADP) National Trends Data (NADP Program Office, 2021). The daily frequency records in SO 2− 4 from these sites were then clipped to match the record length of available DOC and Q data from the associated catchments.

Flow-Adjusted Stream Concentration Data
We flow-adjusted the time series of stream DOC concentration data for each selected catchment to remove the typically dominant Q control from temporal trends and to better isolate potential signals of deterministic trends in stream response (such as stream response to reduced acid deposition). This method has been used in a number of studies that sought to describe temporal trends in stream chemistry constituents that are often confounded by varying stream flow conditions (Hirsch and Slack, 1984;Helsel and Hirsch, 2002;Bekele and Mcfarland, 2004). We first regressed log-transformed DOC concentrations on log-Q data, applying a Locally Weighted Scatterplot Smoothing (LOWESS) algorithm and extracted the residuals. For each regression fit, we used a smoothing pattern coefficient of 0.67 (f = 0.67). One of the benefits of this technique is that a LOWESS model does not assume linearity or normality in the data, therefore not limiting the concentration (C)-Q relationship in each catchment to one specific model (see Supplementary Figure 1 for C-Q plots). We then reordered residuals from the LOWESS fit (i.e., DOC FA concentrations), based on the date/time stamp associated with each observation, and performed a monotonic trend analysis on these residuals using a seasonal Mann-Kendall test (section Seasonal Kendall, see Supplementary Figure 2 and accompanying description for details).

Change Point Detection
We applied a Pettitt test to DOC FA concentrations to identify any threshold effects in trend that would otherwise be missed by simply identifying monotonic trends. The Pettitt test is a nonparametric technique used to identify a single change point in continuous time series data (Pettitt, 1979). This test has been used in a number of studies to detect abrupt changes in hydrologic and climatic variables such as precipitation (Busuioc and Storch, 1996), flood peaks (Liu et al., 2012), and temperature (Wijngaard et al., 2003).  Table 1.

Seasonal Kendall
With both DOC FA and SO 2− 4 deposition data for each catchment, we utilized seasonal Kendall tests to determine the directionality of the aquatic and atmospheric trends. A seasonal Kendall test is widely used as a non-parametric, statistical technique for detecting monotonic trends in time series (Hirsch and Slack, 1984). In this setting, it indicates whether DOC FA , independent of seasonal variability, is increasing or decreasing over time. The product is the normalized test statistic, Kendall's tau, which ranges from −1 to 1. For this analysis, we define "positive trends" as those with a tau value >0.05 and significant beyond an associated alpha threshold of 0.05. "Negative trends" are defined as having a tau value <-0.05 with an alpha threshold of 0.05; and "no trends" are defined as any results that do not meet the requirements of positive or negative trends.

Field Site Descriptions
To test the impact of soil solution and landscape position on DOC dynamics, we sampled soils from two forested headwater catchments with contrasting long-term DOC dynamics for experimentation (Figure 3). The SSHCZO is a 0.08-km 2 forested catchment located in central Pennsylvania in close proximity to Young Woman's Creek (YWC, Figure 3). The catchment is underlain by thick shale and is a first order basin with steep south-and north-facing slopes (25-35%) and narrow ridges. The SSHCZO has a mean annual temperature of ∼9.8 • C and a humid continental climate. Mean annual precipitation averages 1,029-mm and is characterized as acidic (pH ∼4) and stream water is mildly acidic with an average pH of 5.6 (Supplementary Figure 3) (Wen et al., 2020). The soils are derived from colluvium and/or residuum from the underlying shale formation and are Inceptisols and Ultisols with primarily silt loam shallow horizons and a 3 to 5-cm organic layer containing decaying leaf litter (Giardino and Houser, 2015). A number of studies that focused on the cycling of DOC have used the SSHCZO including Andrews et al. (2011), who identified hot spots and hot moments of DOC transport, and Bao et al. (2017), who developed a reactive transport model to predict DOC flux (Andrews et al., 2011;Bao et al., 2017;Wen et al., 2020).
The SRRW, as defined in this study, is a 0.41-km 2 completely forested headwater subcatchment located in eastern Vermont (Figure 3). It is underlain by calcareous granulite and quartz mica phyllite with a forest cover of primarily northern hardwoods such as beech, maple, birch, and ash (Shanley et al., 2004). The catchment encompasses three tributary streams that begin at headwater swamps and move through steep terrain and midelevation benches before meeting in an intersection of slopes at the base of the catchment. The main soil types are Spodosols and Inceptisols in the uplands and Histosols in the lowlands (Kendall et al., 1999). Due to the presence of carbonates in the parent material, ground water and stream water at SRRW are well-buffered with a pH above 7 (Hornbeck et al., 1997;Armfield et al., 2019), also see Supplementary Figure 3.

Soil Core Sampling
In both the SSHCZO and SRRW catchments, soil cores, and samples were collected from the top 10-cm of the soil profile (O/A Horizons) in both swales and planar hillslopes to capture important ecosystem control points. We acquired soil cores for the stop-flow experiment by hammering 2-inch inner diameter PVC pipes with beveled bottoms 10-cm into the soil profile after removing any leaf litter. We then gently dug cores out of the soil profile and capped them to not lose any loose material before being put into a cooler with ice for transport. A total of 18 soil cores were collected for each catchment and sampling time (12 from swales and 6 from planar hillslope locations).
Fall soil samples collected for micro-aggregate separation were acquired from the top 10-cm of the soil horizon in plastic bags and also stored in the cooler for transport before being air-dried. For those experiments designated as occurring in the spring, soil cores were collected from SSHCZO on May 28th, 2019 and SRRW on June 10th, 2019. For the fall experiments, soil cores and samples were collected from SSHCZO on November 8th and SRRW on November 17th.

Experimental Approach to Examining Micro-Scale DOC Dynamics
We conducted stop-flow experiments on soil cores sampled during both the spring and fall seasons to account for seasonal variability within 24-h of sample collection (see Supplementary Figure 4 for a schematic of the core design). Each core was treated with 120-ml of experimental solution that was poured carefully from the top and allowed to saturate the core for 5 min. We used two experimental solutions: the "acidification" solution, simulating acid deposition conditions, had a pH of 3 and an IS of 0.01 M. The "recovery" solution, simulating absence of anthropogenic acids, had a pH of 5.3 and low IS (nanopure water). Note that we refer to this solution as "recovery" solution to indicate near neutral pH and low IS typical for rain in absence of acid deposition and not the ecosystem concept of recovery. After 5 min of interaction time between solution and soils, a valve at the bottom of core was opened, allowing any effluent to drain into a collection vessel. Once the cores drained gravitationally (∼4 min) the valve was closed, and the process repeated three more times yielding 4 separate samples for repeat core treatments. Within 24-h of the experiment, the effluent in each collection vessel was filtered through 0.45-µm polyethersulfone filter before being transferred into combusted amber glass bottles for analysis.
We performed an additional batch experiment with isolated soil micro-aggregates to test specific hypotheses on aggregate break up. To isolate the micro-aggregates, air-dried soil samples were sieved and the fraction smaller than 250-µm and larger than 63-µm was collected. Once isolated, soil aggregates were saturated with one of the two laboratory solutions at a 1:5 aggregate to solution ratio by mass and slowly shaken for 9 min. The aggregates and solution were then gently separated by filtration using a 0.45-µm polyethersulfone filter and a low vacuum pressure.

Laboratory Analyses
To characterize the solid carbon content in SOM for sites without previous data, we analyzed soil subsets of 25 replicated samples collected from SRRW for total organic carbon content (TOC % w/w; note SSHCZO TOC data is available in Andrews et al., 2011). For this, the samples were air-dried, sieved through a 2mm mesh, homogenized using a ball mill and analyzed using a combustion-based elemental analyzer (CE Instruments NC2500) in the Geology Stable Isotope Laboratory at the University of Vermont. Resulting percent carbon values were compared to standards (B2150 for high organic content sediment standard and B2152 for low organic content soil standard) provided by Elemental Microanalysis Limited.
We characterized aggregates from our batch experiments for changes in size, morphology and composition using a VEGA3 TESCAN Scanning Electron Microscope (SEM) and AZtec Elemental Mapping software. We mounted dried aggregates carefully with double-sided carbon tape on metal stubs and sputter-coated them with carbon before visualizing. Because backscattered electron (BSE) SEM mode captures electron density differences well, we used this mode to identify denser mineral grains vs. less dense organic material. BSE was acquired at 5-keV acceleration voltage; and energy-dispersive spectroscopy (EDS) maps were acquired for 5 min with a probe resolution of 15-mm. To compare aggregate size between samples and treatments, we took images at the same magnification (100x) and used image analysis software (ImageJ) to quantify particle size distribution reported in % covered area. This was achieved by converting 100x magnification images of aggregate clusters to binary, scaling, and using the "analyze particles" function with a minimum area set at 1,500 µm 2 and excluding particles that touch the image edges.
We analyzed effluent samples from full soil core experiments for DOC using a Total Organic Carbon Analyzer (Shimadzu, Columbia, MD, USA) within 24-h after collection. To allow for assessment of total DOC release over the entire experiment, we calculated cumulative DOC release; and to allow for comparison between cores, we normalized to solution and soil amount as shown in Equation (1).
To determine significant differences in the cumulative amount of DOC released by categorical independent variables such as treatment and landscape position, Kruskal-Wallis tests were performed. For categorical variables that showed significant results as defined by an alpha threshold of 0.05, a post-hoc Dunn test was used to identify which groups differed from each other group. These statistical tests were also used to determine significant differences in aggregate sizes in the SEM images after the ImageJ particle size distribution analysis was utilized.

Regional Long-Term Trend Analysis
In the Northeastern United States, we identified only nine catchments as having sufficient decadal records of DOC and instantaneous Q to complete our analyses (Table 1 and Figure 4). Of these catchments, six had records that spanned from the 1990s to late 2010s, while three had more limited records between 2003 and 2018 ( Figure 5) Figure 6A). For the time after 2004, data for three additional catchments were available (YWC, CR and WR) and all exhibited negative SO 2− 4 trends ( Figure 6B). Several catchments exhibited increases in DOC FA trends, including catchments experiencing different magnitudes of decrease in wet SO 2− 4 deposition ( Figure 6B). To test for connections between catchment attributes and directionality of DOC FA trends, we plotted Kendall tau values after 2004 against climatic, topographic, and hydrologic attributes from the CAMELS data set (Figure 7). Catchments that have a low number of high precipitation days (i.e., events that are five times greater than the mean daily precipitation) and high % of snow cover, tended to show negative DOC FA trends (Figures 7A,B). In turn, positive DOC FA trends were associated with more days with high precipitation and lower % snow cover. Soil depth varied greatly between catchments and no direct link between DOC FA trends and soil depth was determined ( Figure 7C). However, the link between mean slope and DOC FA trends was more systematic in that catchments with shallow slopes showed no or negative DOC FA trends and catchments with steep slopes showed generally decreasing DOC FA (Figure 7D). Catchments with a low number of days where Q is high (high Q frequency) generally showed positive DOC FA trends while catchments with larger high Q frequency showed negative DOC FA trends ( Figure 7E).

Catchment-Scale Soil Core Experiments at SSHCZO and SRRW
Catchment-scale results of the TOC analysis indicate that carbon content in SRRW soils is variable across landscape positions (Supplementary Figure 5). In the top 10 cm of the swale, TOC was ∼20% as opposed to the planar hillslope, which exhibited ∼5% TOC. The DOC effluent per kg of soil released from soil cores cumulated over the four leaching events was variable across catchments, landscape positions, seasons and treatment. Spring soil cores infiltrated with the recovery-simulating solution released more DOC (14.07 ± 14.27 mg/kg) in comparison to those infiltrated with the acidification solutions (6.95 ± 7.53 mg/kg). Planar hillslopes consistently leached significantly more DOC (17.66 ± 8.57 mg/kg) than swales (4.94 ± 8.98 mg/kg, Figure 8) during this time. For the fall soils, cores that underwent the recovery treatment released more DOC (18.21 ± 15.45 mg/kg) than the acidification (14.46 ± 18.07 mg/kg) treatments and again planar hillslopes released more DOC (17.96 ± 15.10 mg/kg) than the swales (11.89 ± 14.19 mg/kg).

Batch Experiment on Separated Aggregates
The results of the SEM analysis showed differences in aggregate size distribution and morphology by location, landscape position and treatment. The SRRW aggregates in both landscape positions showed heterogeneous associations of larger angular fragments (likely primary minerals) and finer grained materials (Figure 9). Compared to aggregates after acidification treatment, aggregates after recovery treatment appear smaller (e.g., Figures 9B,E). Assessment of particle size distribution from image analyses confirms the qualitative assessment for both swale and hillslope locations. While some larger aggregate/particles are visible after each treatment, the bulk of swale aggregates treated with acidification solution are generally larger (maxima at around 325-µm) than those treated with recovery solution where maximum covered area is at around 175-µm (Figures 9C,F). Aggregate constituents from SSHCZO are homogeneously finegrained and present rounded and tubular particle shapes in both landscape positions (Figure 10). Fine-grained materials (<10µm) are visible for most samples and likely derive from surface coatings and/or disintegrating aggregates. SSHCZO aggregate morphology is similar for both landscape positions and the entire range of particle sizes is represented. Swale aggregates treated with acidification solution show a multimodal particle size distribution with maxima at 225-and 400-µm diameter. After recovery treatment, particles in smaller sizes (below 100µm) are more abundant (Figure 10C). Hillslope particle size distribution is variable and spans the particle size range at similar coverage for both treatments (Figure 10F).

DISCUSSION
Understanding and predicting catchment responses to a regional disturbance is difficult because catchments are spatially heterogeneous systems and exhibit unique characteristics that may engender variable degrees of resistance or resilience to disturbance. This "uniqueness of place" has been acknowledged for many decades in hydrological modeling (Beven, 2000) and the conceptual and numerical integration of single site, single plot observations, stream water patterns and a larger scale driver remains a challenge (Levin, 1992;NSF, 2018). To overcome this uniqueness of place challenge, we suggest an investigative approach that does not seek direct links across scales but that instead evaluates larger-scale pattern and site-specific process in an iterative fashion, and we use reduced acid deposition as test case of regional disturbance for one such iteration. In this context, reduction in acidification in both wet and dry deposition has been hypothesized to have caused an increase in DOC concentrations from many Northern Hemisphere forested catchments (Driscoll et al., 2003;Burns et al., 2006;De Wit et al., 2007;Monteith et al., 2007;Hruška et al., 2009;Sanclements et al., 2018). However, trends in DOC partially contradict each other, which might be due to other and/or overlapping drivers, differences in record length, variable response of unique catchments and spatiotemporal variability. The overarching goal of this study was thus to investigate opportunities and challenges with the process and pattern approach using changes in acid deposition in the Northeastern U.S. as test case.

Pattern-to-Process Investigation: Exploring Catchment Attributes as Covariates in Stream DOC Response to Changes in Atmospheric Deposition
The integration of regional and catchment-scale atmospheric deposition data and catchment attributes (Addor et al., 2017), with co-located stream Q and stream chemistry time-series data provides novel opportunities to investigate catchment response to a regional driver(s). Thus, our first objective was to use our integrated dataset to explore patterns at the regional scale and investigate linkages between reduced acid deposition and stream FIGURE 9 | Secondary electron scanning electron microscopy (SE SEM) images of soil aggregates taken from a swale (A,B) and planar hillslope (D,E) in Sleepers River Research Watershed (SRRW) after treatment with laboratory acidification (A,D) and recovery (B,E) solutions. Particle size distribution from image analyses for swale (C) and planar hillslope (F) aggregates are expressed in % covered area (i.e., the area that is covered by particles that belong to a given size bin).
DOC response and the potential role of catchment characteristics as covariates in that response (objective 1).
To test the connection between decreased atmospheric deposition as a potential driver for DOC stream patterns we compared the directionality of stream water DOC trends against that of SO 2− 4 deposition. The pattern prior to 2004 signals that catchments experiencing the strongest reduction in SO 2− 4 deposition also show the largest increases in DOC FA (Figure 4). However, our results also show that post 2004, the connection between decreases in SO 2− 4 deposition and increases in DOC FA is not evident for each location. For example, the four catchments that show increases in DOC FA post 2004 (SRRW, FB, CR, MB) show large variability in decadal SO 2− 4 deposition trends. These results signal that effects of shifts in atmospheric deposition on DOC are transient or may be confounded by other deterministic trends, and trend analyses thus require explicit temporal contexts. An example are the Catskill catchments (NR and BB), where Burns et al. (2006) found a significant increase in DOC concentrations between 1991 and 2001. Our analysis confirmed this trend for DOC FA until 2004, however, thereafter these catchments exhibit no trends and negative trends (Figure 5). In addition to shifts in directionality of trends for many locations, the Pettitt change point analysis also revealed high variability in timing of the change points even for catchments in close proximity to each other. For example, the shift in DOC FA for BB occurred in 1997, whereas for NR, which drains BB and a larger area, the shift occurs in 2003. As noted in Table 1, these catchments are dramatically different in size, which may in part contribute to the temporal variability in these responses. Furthermore, all these catchments have unique characteristics that likely modulate their response to this disturbance to varying degrees.
To specifically investigate the possible connection between catchment attributes and variable DOC response we investigated the period where DOC FA trends show high inter catchment variability (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). Of all catchment attributes included in the CAMELS dataset (Addor et al., 2017), DOC FA trends show systematic variations with some topographic, climatic, and hydrologic indices. For example, catchments with positive DOC FA trends tend to be flat (i.e., average mean slope is low) while catchments with opposite DOC FA trends have opposite characteristics (Figure 7D). Catchments with the most negative Kendall tau values for DOC FA each have shallow soils; however, shallow soils also produce some of the highest DOC FA values ( Figure 7C). Another finding is that catchments with positive DOC FA trends tend to have less frequent high Q days ( Figure 7E) yet more frequent high precipitation days ( Figure 7A). The association between increasing DOC FA trends and frequency in high precipitation events is in agreement with previous research that indicates DOC export on an annual basis is predominantly driven by large storm/snow melt events (Raymond and Saiers, 2010;Raymond et al., 2016;Wen et al., 2020). Particle size distribution from image analyses for swale (C) and planar hillslope (F) aggregates are expressed in % covered area (i.e., the area that is covered by particles that belong to a given size bin).
Our flow adjustment removes the direct influence of Q on trends in DOC concentrations, thus the apparent link between DOC FA trends and frequency in high precipitation events emphasizes that precipitation drives transport of DOC to streams across catchments with very different characteristics, irrespective of which process led to DOC mobilization. Lastly, catchments that exhibit highest Q frequency also show most negative DOC FA trends ( Figure 7E), which might seem counterintuitive given the fact that the majority of small headwater streams exhibit increasing DOC concentration with increasing Q due to soil flushing (Evans and Davies, 1998;Perdrial et al., 2014;Chorover et al., 2017;Moatar et al., 2017). However, because the flowadjustment process removes the typically dominant Q control on trends, other controls such as reduced acid deposition are emphasized in the data. In this context, flashy catchments (high Q frequency) might more readily deplete carbon pools (Boyer et al., 1997) and lower water residence time, thus leading to decreased DOC FA . In contrast, catchments with less frequent high Q events might access export-controlled DOC pools, especially as the frequency of high precipitation events increased. Our data set is admittedly small, given the limited availability of sites with sufficiently long data records for DOC (Sterle et al., 2019). While our investigation begins to highlight some potential patterns, more research will be necessary to investigate the multi-variate drivers and catchment specific processes that produce such potential patterns.

Process-to-Pattern Investigation: Soil Processes During Simulated Increased and Decreased Acid Deposition Experiments
Our second objective was to investigate the effect of leaching solution composition on select soil processes, specifically aggregate destabilization, as a mechanism for DOC liberation. For this we selected two locations that are well-studied but exhibit contrasting long-term DOC FA trends: SRRW exhibits increasing DOC FA trends, while SSHCZO, a well-studied catchment in close proximity to YWC, shows decreasing DOC FA trends.
In our simulated soil core flushing experiments on top horizons (O-A), soils from SRRW released more DOC than those from SSHCZO for most locations and seasons (Figure 8). These results are in agreement with the higher soil TOC content in SRRW soils (Andrews et al., 2011) and fit the long-term and continued increase in stream DOC over time (Figures 4, 5). Furthermore, SRRW soils show consistently higher DOC releases into recovery solutions, which is again in agreement with the positive DOC FA trend and makes an important connection to reduced acid deposition. The significantly (up to 50%) smaller aggregates after recovery treatment vs. acidification treatment (Figure 9) further corroborate these findings and validate the hypothesized role of aggregates in DOC mobilization during reduced acid deposition in this case. Indeed, previous research on SRRW soils has emphasized DOC release due to the breakup of soil aggregates in lower charge density solutions typical for reduced acid deposition . In previous work the authors note that especially the presence of Ca in SRRW soils might aid aggregation through cation bridging, a process that can be reversed when dilute solutions interact with these aggregates (Armfield et al., 2019;Cincotta et al., 2019).
In contrast, for SSHCZO soil cores, solution chemistry had little or the opposite effect on DOC release (Figure 8), indicating that recovery conditions did not promote DOC mobilization from these soils. Our SEM analyses furthermore indicated little reduction in aggregate size irrespective of the treatment, indicating generally higher aggregate stability (Figure 10). Constituents in SSHCZO aggregates are homogeneously finegrained, which might promote stability through cohesion. The abundant clay mineral content [45%, mostly illite with small amounts of chlorite, vermiculite, and kaolinite (Macdonald et al., 2010)] promotes aggregation (Six et al., 2002) and aids sorption of organic compounds (Kahle et al., 2003). This high carbon stabilization potential might, however, be counteracted through the high water holding capacity that produces favorable condition for microbial processing, leading to overall lower TOC content (Andrews et al., 2011).
For both SRRW and SSHCZO, DOC release varied by landscape position and seasons, irrespective of treatment, and emphasizes the typical spatiotemporal variability in the catchment. For example, planar hillslopes might act as export control points where labile material accumulates until they are hydrologically connected and episodically flushed (Andrews et al., 2011). Our experiments simulate such flushing or export events and mobilize most DOC from planar hillslopes in most cases (Figure 8), despite the fact that TOC content is higher in swales (Andrews et al., 2011). In contrast, swales might represent permanent and/or activated control points where biogeochemical activity is high (but potentially variable) and connectivity to the stream is sustained. Despite ample DOC production, such conditions also lead to losses via microbial processing and/or flushing, thus preventing accumulation (Bernhardt et al., 2017). As a result, cores from swale locations generally release less DOC (Figure 8).
The temporal variability of DOC supply due to seasonal shifts in carbon accumulation, processing and transport is superimposed on the spatial control. For example, both catchments experience spring snowmelt, which flushes soils from otherwise unconnected catchment locations (Boyer et al., 1997) and could lead to a temporal depletion of labile carbon. This effect is visible for SSHCZO, however, the fact that spring SRRW soils release some of the highest amounts of DOC is somewhat counterintuitive: SRRW accumulates a significant snowpack and has a high potential for hydrologic flushing of even distal catchment locations (Shanley et al., 2004;Sebestyen et al., 2008;Armfield et al., 2019). However, it is possible that snowmelt flushing was not extensive enough to deplete stores or that enough time had passed to allow for the replenishing of stocks. For SSHCZO in contrast, we observed the highest DOC release for fall soils, which could be due to the recent addition of leaf litter. Even though leaf litter accumulates readily in swales, these locations again also lose carbon readily (Wen et al., 2020) and the high DOC release from planar hillslopes in fall might therefore again reflect the complex balance between accumulation vs. removal in a catchment setting.

Opportunities and Limitations for the Process and Pattern Investigative Approach
Our experiments emphasize the complexity of catchment specific processes that vary across time (e.g., seasonal dynamics), space (e.g., spatial heterogeneity), or both (e.g., temporal variations in spatial catchment connectivity). Indeed, the differences in DOC release for the SRRW vs. SSHCZO experiments agree with the trend directionality from our regional pattern analyses: tested SRRW soils responded strongly to solutions simulating decreased acid deposition by releasing more DOC and, at a larger scale, SRRW is one of the catchments that shows significant increases in DOC FA . For SSHCZO and the larger trend at YWC, the opposite is the case. In this specific case, observations are consistent across scales, which differs from a similarly designed study on DOC release where lab experiments, field observations and historical data analyses yielded conflicting results (Preston et al., 2011). However, given the inherent complexity of catchment processes, it is unlikely that one experimentally investigated process (in our case aggregate destabilization) explains the variation in long-term trends for catchments in general. To the contrary, this process might be relevant only for catchments of specific combination of biogeophysical and hydrometeorological attributes and results from this work thus provides suggestions for the next iteration of pattern analyses.
For example, wetlands have been reported to play a significant role in DOC dynamics by impacting both water and carbon storage and thus the relative importance (and synchronization) of biogeochemical vs. hydrological processes (Eimers et al., 2008b;Bernhardt et al., 2017;Duan et al., 2017;Kang et al., 2018;Casson et al., 2019). However, while the CAMELS-Chem data set includes soil characteristics that often correlate with wetland cover (e.g., soil conductivity, water fraction, and organic fraction), wetland coverage is not included (Addor et al., 2017) but is a likely driver for our observed differences between sites. Indeed, at the SRRW sites, riparian areas act at least seasonally as wetlands, while at SSHCZO, swale locations are better drained (Andrews et al., 2011;Shanley et al., 2015). Therefore, this work offers opportunities for hypothesis generation and further process investigations for catchments with contrasting longterm patterns.
Integrated datasets, such as the CAMELS-Chem dataset used in this study, are fundamental to this type of research, and continued and concerted efforts to generate and compile these types of data in readily-accessible formats will help address current limitations for cross-scale investigations. For example, only ∼7% of CAMELS catchments had DOC records sufficiently long and consistent enough to perform robust long-term trend analyses, with only ∼2% located in the Northeastern US. The CAMELS data set includes a wide variety of catchment attributes; however, it is likely that important attributes are not captured. These spatial and temporal gaps currently limit statistical power and the development of a larger consensus on processes, suggesting a need for more long-term monitoring and data fusion/integration efforts as illustrated here. Overall, consistent, extensive, and multiple-perspective monitoring systems are urgently needed to record long-term alteration of water quantity and quality response to climate change and human perturbation efforts (Lovett et al., 2007;Magner and Brooks, 2008;Li et al., 2021;Zhi et al., 2021). In the meantime, we believe that such integrated datasets, novel data science tools, and process investigations will allow the catchment science community to make progress addressing the problem of scale.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: the dataset [CAMELS-Chem] for this study will be made available in Sterle et al. (2019) upon publication. The stream water dataset analyzed for DOC trends at SRRW can be found at: doi: 10.5066/P9380HQG and doi: 10.5066/P929KMVK.

AUTHOR CONTRIBUTIONS
TA was responsible for sample collection, preparation, analyses, and production of both figures and tables. JP aided in interpretation, authorship, field work, lab work, lab training, and allowed lab use. KU, DR, AH, AL, and LL provided essential expertise to the interpretation of the results. GS compiled and collected data used in analyses. LS, CB, and HW aided in field work, laboratory work, and critical analyses. NP contributed through interpretations, authorship, and lab training. All authors contributed to the article and approved the submitted version.