Flow Partitioning Modelling Using High-Resolution Isotopic and Electrical Conductivity Data

Abstract

Water-stable isotopic (WSI) data are widely used in hydrological modelling investigations. However, the long-term monitoring of these tracers at high-temporal resolution (sub-hourly) remains challenging due to technical and financial limitations. Thus, alternative tracers that allow continuous high-frequency monitoring for identifying fast-occurring hydrological processes via numerical simulations are needed. We used a flexible numerical flow-partitioning model (TraSPAN) that simulates tracer mass balance and water flux response to investigate the relative contributions of event (new) and pre-event (old) water fractions to total runoff. We tested four TraSPAN structures that represent different hydrological functioning to simulate storm flow partitioning for an event in a headwater forested temperate catchment in Western, Oregon, USA using four-hour WSI and 0.25-h electrical conductivity (EC) data. Our results showed strong fits of the water flux and tracer signals and a remarkable level of agreement of flow partitioning proportions and overall process-based hydrological understanding when the model was calibrated using either tracer. In both cases, the best model of the rainstorm event indicated that the proportion of effective precipitation routed as event water varies over time and that water is stored and routed through two reservoir pairs for event and pre-event. Our results provide great promise for the use of sub-hourly monitored EC as an alternative tracer to WSI in hydrological modelling applications that require long-term high-resolution data to investigate non-stationarities in hydrological systems.

1. Introduction

The identification of the water sources contributing to runoff is fundamental to understand the linkage and interactions between water and biogeochemical cycles and the transport of contaminants and solutes at the catchment and landscape scales [1,2,3,4]. During precipitation events total runoff can be partitioned into event—“new” water from incoming precipitation—and pre-event water—“old” water stored in the catchment prior to a given precipitation event [5,6]. Depending on the catchment conditions (e.g., vegetation, soil type, geology, topography, antecedent moisture) and event characteristics (e.g., precipitation amount and temporal variability) the event and pre-event water fractions vary. As such, understanding the spatial and temporal variability of the contributions of different water pools to the hydrograph is not only a fundamental question in hydrological science [7], but is also needed for the implementation of effective water resources management strategies worldwide [8].

Given that the contribution of different water pools to total runoff is time dependent owing to the time variant nature of flow response processes [7,9,10,11,12], we struggle to apply effective monitoring strategies in catchments with different environmental conditions [13]. Conservative (e.g., water stable isotopes, ²H and ¹⁸O, and chloride) and non-conservative (e.g., electrical conductivity (EC) and silica) tracers have been used to constrain mixing models of flow partitioning at the event scale (e.g., [14,15,16,17,18]). One of the most commonly applied methods is tracer-based two-component end member mixing analysis [5,6,19,20] to quantify the proportions of event and pre-event water to total runoff using a simple mass balance (cf., Buttle [5] and Klaus and Mcdonnell [6] for reviews). However, since these models only account for the mixing of the tracer within a hydrological system, they provide only limited information about the combined flow and tracer mixing dynamics in response to precipitation inputs. Therefore, these models have limited ability to provide a process-based understanding of catchment behavior [21,22,23].

In response to this challenge, numerical models in which the tracer mass balance and the hydrological flow response are coupled have been developed [24,25,26,27]. These models allow for the simultaneous simulation of the streamflow hydrograph—water flux—and the tracer mixing. In general, these models account for the water flux partitioning by applying the unit hydrograph approach, and the tracer mixing using travel time distributions (TTDs) [28]. The application of TTDs accounts for the estimation of the possible travel times of the tracer within the system. That is, the time that a water molecule takes to travel within a hydrological system from the moment it enters as precipitation or snowmelt to the time it exists as runoff [29,30]. The TTDs shapes are used to identify runoff processes [7,29,31,32] by providing information about the physical processes that influence the internal mixing of different water sources.

Despite the advantages of numerical approaches for hydrograph separation, the availability of long-term high-resolution water geochemical data remains a challenge. To date, the majority of applications of tracer-based hydrograph separation techniques have been conducted using WSIs as tracers [6]. Despite their recognized usefulness and reliability (being conservative), WSIs sampling at fine resolution (e.g., sub-hourly) is still sparse given high associated costs and, thus, limits the description of the rapid response of streamflow to water inputs and the inter-storm variation of the input isotopic composition [33,34,35,36]. This limitation has resulted in high uncertainties in the estimation of flow components [35,37]. Recent studies have monitored WSIs at high-temporal resolution (every 30 min) [38,39,40,41,42,43] and applied simple mixing models to partition flow components [9,18,44,45]. However, while technological developments currently allow the deployment of field analyzers to measure the WSI of inputs and outputs (e.g., rainfall and streamflow) at high temporal resolution it is unfeasible to broadly implement such analyzers. Thus, there is a need for alternative inexpensive and low maintenance water quality parameters (i.e., tracers) that allow investigating internal catchment processes at a high resolution [13,46].

Electrical conductivity, or specific conductance (EC) of water, is an alternative tracer often used in flow partitioning and water quality studies, either alone or in combination with WSIs (e.g., [16,20,47,48,49,50,51,52,53,54]). The main advantage of using EC is that it can be continuously monitored at high temporal resolution (seconds to minutes) using inexpensive in-stream probes [13,46,48]. Despite the non-conservative nature of EC as it highly depends on the water contact time with the mineral substrate, in particular [53,55,56], EC has yielded similar streamflow portioning results than WSIs using traditional mass balance models [16,51]. However, its effectiveness has not been tested against WSI using sophisticated hydrological modeling approaches.

In this study, we compare the results of a flow partitioning tracer-based hydrological model calibrated using WSIs and EC. Our specific objectives are: (1) to evaluate different model structures (representing different assumptions of the internal catchment hydrological functioning) and determine the model structure that best simulates flow partitioning using WSI and EC data; and (2) to compare the process-based hydrological understanding obtained from the models calibrated using each tracer.

2. Materials and Methods

2.1. Study Area

This study was conducted at the Mack Creek catchment (5.8 km²) in the H.J. Andrews Experimental Forest in the Western Cascades of Oregon (Figure 1). Mack Creek is a tributary of Lookout Creek, which drains into the Blue and McKenzie rivers within the Willamette River Basin. Glaciation occurred in the catchment leaving U-shaped valley morphologies with a steep slope (average 46%). The elevation ranges between 758 and 1610 m above sea level (a.s.l.) [57,58]. The catchment is underlie by ridge-capping andesite lava flows (upper Sardine Formation and Pliocene flows). The soils in the catchment are gravelly loams at high elevation that transition into the predominantly gravelly sandy loams that characterize over 70% of Mack Creek drainage area [59]. The forest is dominated by 400–500 year-old coniferous trees, including Douglas fir (Psuedotsuga menziesii), western hemlock (Tsuga heterophylla), and western red cedar (Thuja plicata). The climate is Mediterranean with wet, mild winters and dry summers. Fall and winter precipitation falls as a mix of rain and snow and may accumulate and last from early November to late June [60]. Mean annual precipitation (between 2002 and 2016) was 2243 and 2709 mm at the CENMET (1018 m a.s.l.) and UPMET (1294 m a.s.l.) meteorological stations (Figure 1), respectively. However, precipitation in 2015—when we conducted this study—was only 51% of this 15-year average and fell almost entirely in the form of rain.

2.2. Hydrometric Data Collection

We monitored water fluxes—precipitation (P) and streamflow (Q)—during a large rainstorm event between 28 October and 7 November 2015. We used 5-min resolution P and Q data maintained as part of the NSF Long-Term Ecological Research (LTER) program [61]. We used Q data derived from an 18-in flume (GSMACK; Figure 1) [61] and mean P values across the UPMET station and a precipitation tipping bucket located in the roof top of the GSMACK [62]. The 8-in diameter GSMACK gauge is heated with a propane heater. Continuous precipitation is recorded with depth measurements using a Stevens A-35 chart recorder or a Stevens PAT water level shaft encoder. The UPMET precipitation gauge is stand-alone composed of standing pipe with tank gauge, a propane-heated 20-inch diameter orifice, surrounded by a Valdai-style double wind fence. A temperature probe controls the orifice heating by turning a pump and heater on/off. The stand-alone rain gauge was specifically developed to withstand heavy snow with depths up to 3–4 m and windy condition.

2.3. Water Stable Isotope Data Collection and Analysis

We collected 44,500 mL streamflow grab samples every 4 h using an automatic autosampler (ISCO-3700). A grab sample was also collected before the event in order to characterize the streamflow isotopic composition during base flow conditions (i.e., pre-event). We described the event isotopic signature of precipitation at a meteorological station located at 430 m a.s.l. (PRIMET; Figure 1) with a sequential rainfall sampler designed to collect up to 12,200-mL water samples to characterize the temporal variability of water stable isotopes [33,63]. The sequential collector was located in a clearing, thus, we assumed that the difference between the isotopic signature of direct rainfall and throughfall was minimal [64]. All water samples for isotopic analysis were collected and stored in 20-mL glass bottles with conical inserts without headspace and kept in dark at relatively cool conditions (~15 °C) to avoid fractionation by evaporation until their analysis was conducted at the Watershed Processes Laboratory at Oregon State University.

We measured the water-stable isotopes (δ¹⁸O) in all samples using a cavity ring down spectroscopy liquid and vapor isotopic analyzer (Picarro L2130-i, Picarro Inc, Santa Clara, CA). We ran all samples under high-precision, including six injections per sample. The first three injections were discarded to account for memory effects. Two internal (secondary) standards (MET1, δ¹⁸O = −14.49‰ and BB1, δ¹⁸O = −7.61‰) were used to develop calibration curves, while a third internal standard (ALASKA1, δ¹⁸O = −11.09‰) was used to estimate a drift correction equation (drift was always below 0.000152‰). All internal standards were calibrated against the IAEA primary standards for the Vienna Standard Mean Ocean Water (VSMOW2, δ¹⁸O = 0.0‰), Greenland Ice Sheet Precipitation (GISP, δ¹⁸O = −24.76‰), and Standard Light Antarctic Precipitation (SLAP2, δ¹⁸O = −55.5‰). The uncertainty in our secondary standards (i.e., standard deviation) is <0.01‰. Based on >50 duplicate samples (collected concurrently under comparable conditions) in rainfall and streamflow, we estimated an internal laboratory precision of 0.03‰. The external accuracy of our laboratory was 0.06‰. This accuracy was computed as the mean difference between 60 estimated values and a known water standard. We used only the internal precision as an overall measure of accuracy although we acknowledge that there are other sources of uncertainty [65].

2.4. Electrical Conductivity Measurements

We measured the conductance of the sequential water samples collected for the water stable isotopes analysis using a Hanna^® Multiparameter (HI9828) Water Quality Portable Meter and used these measurements to characterize the EC of the precipitation during the event. The portable meter was calibrated in the laboratory following manufacture guidelines. The manufacture accuracy of this instrument is 1% or 1 μS/cm (whichever is larger). Replicate measurements for 30 samples indicated a precision of <5%. The EC measurements were conducted after the 20-mL samples for water stable isotopes were collected to avoid sample contamination. Streamflow EC was continuously recorded (every 5 min) at the GSMACK gauge with a Campbell Scientific CS547A EC and temperature probe (with a 5% accuracy) [66]. The EC data is corrected with a temperature coefficient based on YSI Pro30 conductivity instrument [66].

2.5. Tracer-Based Hydrograph Separation Modeling

In recent years, different approaches that allow the incorporation of geochemical tracers mixing into hydrological models have been developed. The application of conceptual models using TTD functions [25,27,67] and ranked storage-age-selection (rSAS) functions [68] are amongst the most frequently applied approaches. We selected the former approach as it allows for the implementation of different model structures (representing different catchment hydrological behavior) that have been widely tested using WSI for model calibration in other studies and, thus, allow for a direct comparison with our results using both high-resolution WSIs and EC data.

We used the flexible modelling framework developed by Segura et al. [25] that considers both the tracer mass balance and the hydrological flow response during rainfall-runoff events. We refer to this model as TraSPAN (Tracer-based Streamflow Partitioning ANalysis model). TraSPAN shares similarities with other tracer-based hydrological models (e.g., [24,27]). It allows for the implementation of different model structures to simulate different levels of hydrological complexity, assuming that rainfall-runoff can be partitioned into event and pre-event components. As such, we implemented this framework to evaluate competing hypotheses of hydrological response against observed Q and tracer data to select the best model structure for the given hydrological system [26,69,70].

TRaSPAN is composed of three modules (

Table 1. Modules, equations, and parameters of the four Tracer-Based Streamflow Partitioning Analysis (TraSPAN) model structures.

Module		Parameters and Equations		Units	Model Structure
					1	2	3	4
1 a	Peff	Initial antecedent rainfall index, So		-	X	X	X	X
		Memory timescale parameter, w		time steps d	X	X	X	X
	Equations	$P_{eff}\left(t\right)=p\left(t\right)s\left(t\right)$	(1)
		$s\left(t\right)=c·p\left(t\right)-s\left(t-∆t\right)\left(1-\frac{1}{w}\right)$	(2)
		where p is precipitation and c is the normalization constant to maintain ${\displaystyle \sum }{P}_{eff}={\displaystyle \sum }Q$
2 b	f, constant	Fraction of effective rainfall routed as event water, f		-	X	X
	f, variable	Normalization constant, cf		15-min/mm			X	X
		Memory timescale parameter, wf		time steps d			X	X
	Equations	$f\left(t\right)=c_f·p\left(t\right)-f\left(t-∆t\right)\left(1-\frac{1}{w_f}\right)$	(3)
3 c	Single Reservoir	Mean transit time of event water, ke		h	X		X
		Time delay for the event fraction response, elag		h	X		X
		Mean transit time of pre-event water kp		h	X		X
		Time delay for the pre-event fraction response, plag		h	X		X
	Equations	$Q_e\left(t\right)={{\displaystyle \int }}_0^t{h}_e\left(\tau \right)P_{eff}\left(t-\tau \right)f\left(t-\tau \right)d\tau$	(4)
		$Q_p\left(t\right)={{\displaystyle \int }}_0^t{h}_p\left(\tau \right)P_{eff}\left(t-\tau \right)\left[1-f\left(t-\tau \right)\right]d\tau$	(5)
		$he\left(\tau \right)=e_{lag}-\frac{1}{k_e}\exp \left(\frac{\tau }{k_e}\right)$	(6)
		$hp\left(\tau \right)=p_{lag}-\frac{1}{k_p}\exp \left(\frac{\tau }{k_p}\right)$	(7)
		$C_t\left(t\right)=\frac{Q_e\left(t\right) C_e+Q_p\left(t\right) C_p}{Q_e+Q_b}$	(8)
		where Ce is the tracer composition of event water (rainfall) and Cp is the tracer composition in pre-event water (e.g., baseflow prior to the rainstorm).
	Two Parallel Reservoirs	Fraction of water routed into the fast-responding reservoir of event water, qe		-		X		X
		Mean transit time of the fast fraction of event water, kfe		h		X		X
		Mean transit time of the slow fraction of event water, kse		h		X		X
		Time delay for the event fraction response, elag		h
		Fraction of water routed into the fast-responding reservoir of pre-event water, qp		-		X		X
		Mean transit time of the fast fraction of pre-event water, kfp		h		X		X
		Mean transit time of the slow fraction of pre-event water, kse		h		X		X
		Time delay for the pre-event fraction response, plag		h		X		X
	Equations	$he\left(\tau \right)=-\frac{q_e}{k_{fe}}\exp \left(-\frac{\tau -e_{lag}}{k_{fe}}\right)+\frac{1-q_e}{k_{se}}\exp \left(\frac{\tau -e_{lag}}{k_{se}}\right)$	(9)
		$hp\left(\tau \right)=\frac{q_p}{k_{fp}}\exp \left(-\frac{\tau -p_{lag}}{k_{fp}}\right)+\frac{1-q_p}{k_{sp}}\exp \left(\frac{\tau -p_{lag}}{k_{sp}}\right)$	(10)
		Equations (4), (5) and (8) are also used here.
		Total number of parameters			7	11	8	12

^a Module 1: Effective rainfall module; ^b Module 2: Event and pre-event routing module; ^c Module 3: Event and pre-event transfer functions module; ^d 15 min time steps.

). Module 1: The effective rainfall module is a non-linear routine that computes the effective rainfall (P_eff) as the product of precipitation (P) and the antecedent rainfall index (s) [71] (Equations (1) and (2) in Table 1). Module 2: The event and pre-event routing module defines the fraction P_eff routed as either event or pre-event water. Module 3: The Event and pre-event transfer functions module includes the TTDs to convolve the fractions of event and pre-event water and, thus, represents the internal hydrological behavior of the system. Modules 2 and 3 are flexible, allowing for the fraction (f) of P_eff routed as event water to be constant or time-variant (i.e., to vary over the duration of the event; Equation (3) in Table 1) and for the incorporation of TTDs of varying degrees of complexity into the convolution integrals for event ant pre-event water (Equations (4) and (5)). Even though there exists a variety of theoretical TTDs that could represent tracer mixing in hydrological models, we selected two of the most commonly applied TTDs that have been tested to represent hydrological systems in other catchments [27,72,73]. These are the exponential model (EM; [74]) representing a single linear reservoir (Equations (6) and (7) in Table 1) and the two-parallel linear reservoirs model (TPLR; [27]) representing two connected linear reservoirs (Equations (9) and (10) in Table 1). A detailed description of the model framework and its modules can be found in Segura et al. [25]. We also allowed the tracer signal to vary 10% around the concentration measured before the event to account for the uncertainty in pre-event water composition.

The fractions of event and pre-event water are routed according to the given TTD for h_e(τ) and h_p(τ) (Table 1, Figure 2) and the resulting event (Q_e) and pre-event (Q_p) fractions are used to calculate the tracer concentration based on a mass balance approach (Equation (8), Table 1). Total discharge (Q_t) is the sum of event and pre-event water plus the baseflow (Q_b), which was subtracted from the input data prior to the modelling. We defined Q_b as the discharge prior to the beginning of the rainfall event. All the parameters in the model were estimated by the simultaneous calibration of the streamflow and tracer data [25].

We explored four TraSPAN model structures of varying complexity depending on the treatment of the fraction of P_eff routed as event water and the number of reservoirs used to route the event and pre-event water fractions (Figure 2). The number of parameters fitted in each of the four model structures varied between 7 and 12 (Table 1). Structure 1 had a constant fraction f of P_eff routed as event water and a single reservoir for each of the event and pre-event water components (seven parameters). Structure 2 also had a constant fraction f, but two reservoirs in parallel for routing the event and pre-event water components (11 parameters). Structure 3 had a time-variant fraction f(t) of P_eff routed as event and pre-event water and a single reservoir for each of the event and pre-event water components (8 parameters). Structure 4 also had a time-variant fraction f(t), but two reservoirs in parallel for each of the event and pre-event water components (12 parameters). For all the model structures, Q, δ¹⁸O, and EC data were computed every 15 min.

2.6. Model Simulations, Performance, and Uncertainty

Initially, we ran each of the model structures per tracer 80 million times within a broad range of Monte Carlo-generated parameter values from uniform distributions (

Table 2. Parameter sets that yielded the highest NSEs values for each TraSPAN model structure and tracer used for calibration.

Model Structure	Tracer		Module1		Module 2			Module 3
			So	w	f	cf	wf	ke (h)	elag (h)	kp (h)	Plag (h)	qe	kfe (h)	kse (h)	qp	kfp (h)	ksp (h)
1	δ18O	Min	0.00	1.00	0.17			15.06	0.00	13.98	0.00
		Max	1.00	400.00	0.33			55.84	5.75	74.26	3.75
		Best	0.99	4.03	0.28			28.18	0.25	34.07	0.00
	EC	Min	0.00	11.03	0.16			15.96	0.00	12.67	0.00
		Max	0.75	399.98	0.24			50.53	3.75	49.47	2.75
		Best	0.01	390.44	0.20			21.19	0.25	19.64	0.25
2	δ18O	Min	0.00	1.00	0.15				0.00		0.00	0.00	0.03	16.24	0.00	0.03	6.31
		Max	1.00	400.00	0.37				6.00		4.50	1.00	49.87	250.0	1.00	49.99	250.0
		Best	0.65	6.60	0.25				0.75		0.25	0.66	16.66	128.7	0.15	6.35	47.23
	EC	Min	0.00	1.09	0.13				0.00		0.00	0.00	0.25	17.82	0.00	0.11	5.75
		Max	1.00	400.00	0.27				5.25		4.25	1.00	47.78	250.0	1.00	46.99	250.0
		Best	0.01	190.67	0.19				0.75		1.00	0.69	6.18	122.1	0.66	4.89	208.8
3	δ18O	Min	0.00	1.00		0.01	1.01	17.83	0.00	18.23	0.00
		Max	1.00	39.99		0.19	39.99	43.40	2.00	72.50	3.75
		Best	0.08	24.90		0.02	21.33	26.10	0.00	30.61	0.00
	EC	Min	0.00	1.00		0.01	1.00	22.80	0.00	17.50	0.00
		Max	1.00	40.00		0.15	40.00	64.50	6.50	64.84	3.25
		Best	0.75	2.50		0.02	12.45	37.63	0.25	35.53	0.25
4	δ18O	Min	0.00	1.02		0.01	1.23		0.00		0.00	0.06	1.06	22.9	0.02	1.25	20.77
		Max	1.00	39.98		0.14	39.99		2.25		3.25	1.00	31.30	249.8	1.00	49.54	250.0
		Best	0.31	39.63		0.01	34.95		1.00		0.50	0.48	6.24	45.00	0.31	6.57	44.26
	EC	Min	0.00	1.01		0.01	1.00		0.00		0.00	0.00	0.24	27.2	0.00	0.04	13.60
		Max	1.00	40.00		0.16	39.98		3.25		3.00	1.00	49.95	250.0	1.00	38.49	249.9
		Best	0.11	21.03		0.01	32.75		0.00		0.00	0.63	20.62	140.6	0.27	3.84	59.94

). We used the Nash-Sutcliffe efficiency (NSE) coefficient [75] during the calibration always considering the average NSE for the fits of Q and either δ¹⁸O or EC. In order to improve the identification of parameters we ran each model structure 20 million additional times using a narrower range of parameters corresponding to those that yielded an NSE of at least 80% of the maximum NSE during the first run. Given that structure 4 had a larger number of calibrated parameters, we conducted a third model run for 10 million additional times considering the values of the parameters that yielded NSE values of at least 90% of the maximum NSE of the second run as initial parameter ranges [76].

Overall model performance was evaluated using the NSE as a measure of goodness of fit, and the Akaike information criterion (AIC) [77]—using the χ² statistic as the likelihood [78]—as a measure of model parsimoniousness. This was necessary considering the different number of parameters of the considered model structures (Table 1). We considered the sum of AICs for the Q and tracer data. Given that the χ² test requires the uncertainty in the observations, we quantified it for Q, δ¹⁸O, and EC. The uncertainty in the discharge was calculated as the sum of the uncertainty associated with the stage measurements and the uncertainty associated with the rating curve. The stage uncertainty was estimated to be less than 1.2 mm considering both the instrumentation precision and observed bias. We assessed the uncertainty in the rating curve by conducting 10,000 Monte Carlo simulations, in which its parameters were randomly varied. This uncertainty varied between 0.9% and 4.31% with an average of 3.2%. The total Q resulting uncertainty for our monitored storm varied between 2% and 14%. The uncertainty in δ¹⁸O was estimated to be 0.03‰ (Section 3.4) and the uncertainty in EC was assumed to be 5%, considering the precision of the instrument.

We evaluated the uncertainty in the model parameters based on the results of their last run using a threshold of behavioral solutions. This threshold of behavioral solutions was set up to include parameter sets that yielded NSEs above 0.70 (for model structures 1 and 3) and above 0.75 (for model structures 2 and 4). In all cases, there were >5000 behavioral sets. The difference in threshold results from the larger number of simulations for models 1 and 3 that would have been required in order to get 5000 parameter sets with NSE > 0.75.

2.7. Comparison of Model Simulations for δ¹⁸O and EC

Once the model structure that best simulated the observed Q and tracer data when calibrated for both tracers (δ¹⁸O and EC) was selected, we compared the results in terms of their capability to simulate the Q and tracer data observations, their distribution and ranges of behavioral parameters, and their TTDs. The TTD comparison, which provides information about the mixing and transport processes within the hydrological system, allowed for a direct comparison of the hydrological behavior determined by the model structures using each tracer.

3. Results

3.1. Hydrometric and Tracer Characterization

The monitored rainfall-runoff event lasted for 10 days and 6.5 h (t₀ − t_f in Figure 3a). During this period, total precipitation (in the form of rain; P) and total runoff (Q) reached 155.1 mm and 20.6 mm, respectively. These values corresponded to a runoff coefficient (Q/P) of 0.13. During the first 44.5 h (t₀ − t₁), total P was 19.9 mm with a mean intensity of 0.45 mm h⁻¹. This period was characterized by little response in the hydrograph. Subsequently (t₁ − t₂), Q started to increase in response to P inputs of higher intensity (0.95 mm h⁻¹). At t₂ (92.5 h since t₀) Q reached its maximum after a total P amount of 107.4 mm. Later (t₂ − t₃), the recession of the hydrograph started as P intensity decreased to 0.30 mm h⁻¹. Then, at t₃, the last Q peak occurred 112 h after the beginning of the event. By this time, total P was 140.7 mm. Between t₃ and the end of the event (t_f), P almost completely ceased, and the hydrograph recession proceeded. During this last period (t₃ − t_f), total P was 14.4 mm with a mean intensity of 0.06 mm h⁻¹.

The temporal dynamics of the δ¹⁸O and EC signals of P and Q are shown in Figure 3b,c respectively. The δ¹⁸O of P (δ¹⁸O_P) at t₀ (baseline signal) was the lowest during the event (−10.66 ‰) and started increasing as P inputs increased (Figure 3b). The δ¹⁸O_P signal peaked at −3.5‰ in the water sample collected a few hours before the Q peak at t₃. After this time, during the recession of the hydrograph, the δ¹⁸O_P decreased to the baseline value in the sample collected at the end of the event (t_f). Similarly, the EC signal in P (EC_P) started at a low baseline value of 5.2 μS cm⁻¹, and increased as P inputs increased (Figure 3c). The EC_P peaked at a value of 30 μS cm⁻¹ in the sample collected a few hours before the highest peak of the hydrograph at t₂. After this period, EC_P decreased near the baseline value, and then peaked again (as δ¹⁸O_P did) in the sample collected before the last Q peak (around t₃) with a value of 11 μS cm⁻¹. Then, EC_P decreased to the baseline value during the recession of the Q as P inputs decreased until the end of the event (t_f). The temporal variability in the Q δ¹⁸O (δ¹⁸O_Q) was smaller (−10.92‰ to −8.3‰; Figure 3b) than the variability observed in δ¹⁸O_P. The δ¹⁸O_Q showed a pattern related to the hydrograph, i.e., δ¹⁸O_Q increased as Q increased, and vice versa, in response to the dynamics of P. Similarly, the temporal change of EC in Q (EC_Q) was lower (31.8–41.6 μS cm⁻¹; Figure 3c) than EC_P. The EC_Q variability was the mirror image of the δ¹⁸O_Q signal, i.e., when the δ¹⁸O_Q signal (and Q) increased, the EC_Q signal decreased, and vice versa (the r² between EC and δ¹⁸O. is 0.69). The Q values of the pre-event water fraction were −10.92‰ for δ¹⁸O and 41.0 μS for EC.

3.2. Model Performance

A summary of the parameter values that yielded the highest NSEs and the statistics of the calibration model performance for the four structures using each tracer are presented in Table 2 and

Table 3. TraSPAN model performance in terms of the mean Nash-Sutcliffe efficiency (NSE) for discharge and tracer concentrations and the Akaike Information Criterion for discharge (AIC_Q) and tracer (AIC_T, i.e., δ¹⁸O or EC) values.

Tracer	Model Structure	Pre-Event Fraction (%)	NSE	AICQ	AICT	ΣAIC *
δ18O	1	71.7	0.86	20,762	644	21,405
	2	76.2	0.87	16,987	514	17,501
	3	76.3	0.81	21,563	1025	22,589
	4	74.9	0.87	11,390	849	12,238
EC	1	79.6	0.79	20,762	448	21,209
	2	79.5	0.86	15,517	379	15,896
	3	81.7	0.86	22,581	213	22,793
	4	81.0	0.90	10,747	245	10,992

Note: Bold values indicate the best model structure considering NSE or AIC. * ΣAIC = AIC_Q + AIC_T.

respectively. In general, all model structures using both tracers yielded strong fits (NSE > 0.79) to the observed Q and tracer data (Table 3). The highest NSEs were obtained using model structures 2 and 4 (NSE = 0.87) for δ¹⁸O and model structure 4 (NSE = 0.90) for EC. Our evaluation of the models structures’ parsimoniousness showed that the variation of the AIC values for Q (10,747–22,581) was higher than the variation of the AIC values for the tracers (213–1025) among the different model structures (Table 3). Even though model structure 4 had the largest number of fitting parameters, the sum of the AIC values for Q and the tracers (ΣAIC) showed that this model structure yielded the lowest values for both δ¹⁸O (ΣAIC = 12,238) and EC (ΣAIC = 10,992).

3.3. Modelled Streamflow Partitioning

The fitted hydrographs using TRaSPAN model structure 1 (Figure 2a) for δ¹⁸O and EC (Figure 4a,d) poorly resembled the temporal dynamics of Q during the monitored event. In general, the model was not able to reproduce the Q peak responses to P inputs, or the recession limbs of the hydrographs. This was particularly noticeable during the highest Q peak and rapid recession of the hydrograph on DOY 305 and by the overestimation (DOY 306–308) and underestimation (DOY 308–311) of Q during the recession limb of the hydrograph. Similarly, this model structure poorly reproduced the δ¹⁸O_Q and EC_Q temporal variability during the event. The structure could not resemble the δ¹⁸O_Q tracer signal dynamics at the beginning of the event (DOY 301–303) and the enriched isotopic composition during the highest Q peak response on DOY 306 (Figure 4b). During the remaining of the event, the tracer signal was better captured within the uncertainty bands of the simulations. Regarding the EC_Q temporal dynamics, the model simulations were even weaker (Figure 4e). The model poorly reproduced this tracer’s signal at the beginning of the event (DOY 301–304), during the highest Q peak response (DOY 305), and during the last part of the event’s recession limb (DOY 311). The pre-event water fractions estimated using this model structure were 71.7% for δ¹⁸O and 79.6% for EC (Table 3) and the temporal dynamics of the event and pre-event water contributions to total Q are shown in Figure 4c,f, respectively. These fractions depicted similar overall model results using both tracers. Our simulations also indicated an over prediction of the pre-event water fraction (i.e., Q_e > 100%) at the beginning of the event with δ¹⁸O (DOY 301–302; Figure 4c) and at the first part of the recession limb (DOY 306–308; Figure 4d) with EC. Similar results were found for model structure 3 (time-variant fraction of P_eff routed as event water and single reservoirs for event and pre-event; Figure 2c). For this structure, the implementation of the time-variant routine for the fraction of P_eff routed as event water did not improve the simulations of the hydrograph and the temporal variation of the tracers (Figure S2). The pre-event water fractions estimated using this model structure were 76.3% for δ¹⁸O and 81.7% for EC (Table 3) and the temporal dynamics of the event and pre-event water contributions to total Q were similar using both tracers (Figure S2c,f).

The simulation results using TRaSPAN model structure 4 (i.e., the time variant fraction of P_eff routed as event water and two connected reservoirs for event and pre-event water; Figure 2d) showed that this structure better captured the temporal dynamics of Q using both tracers. All peaks (e.g., the highest peak produced on DOY 305) and the recession (DOY 306–311) of the hydrograph were well simulated by this model structure using δ¹⁸O (Figure 5a) and EC (Figure 5d). Similarly, the temporal dynamics of the tracers’ signal were well simulated by this model structure. Even though this model structure underestimated the δ¹⁸O_Q enrichment (Figure 5b) and overestimated the decrease in EC_Q (Figure 5e) on DOY 304; the most enriched δ¹⁸O_Q and the lowest EC_Q peak values (DOY 305), as well as the rest of the tracer dynamics, fitted the observed data well. The estimated proportions of pre-event water were 74.9% using δ¹⁸O and 81% using EC (Table 3). The temporal variability of the pre-event water contributions were similar using δ¹⁸O (Figure 5c) and EC (Figure 5f) and showed a dominance of pre-event water contributions, even during the peak Q generation.

Even though structure 2 (i.e., the constant fraction of P_eff routed as event water and two connected reservoirs for event and pre-event; Figure 2b), simulated the hydrograph better than structures 1 and 3 using both tracers (e.g., it captured the Q peaks better), it still had issues simulating the peak hydrograph response compared to structure 4, particularly when calibrated using δ¹⁸O (Figure S1a,b). Structure 2 also had issues replicating the δ¹⁸O_Q (Figure S1c) and EC_Q (Figure S1d), particularly when calibrated using EC during the rising limb of the hydrograph. The pre-event water fractions estimated using this model structure were 76.2% for δ¹⁸O and 79.5% for EC (Table 3). However, in contrast to all of the other model structures, the temporal dynamics of the event and pre-event water contributions to total Q were different using both tracers. This structure tended to estimate lower pre-event water contributions during the rising limb of the hydrograph and higher pre-event water contributions during the recession when calibrated using δ¹⁸O (Figure S1c). An opposite trend was observed when calibrated for EC (Figure S1f).

3.4. Model Parameter Identification for the Best Model Structure

According to the dotty plots (i.e., values of the NSE as a function of the calibrated parameters for a given model run), seven of the 12 parameters of model structure 4 using EC for calibration reached a single peak in the parameter space yielding the highest NSE (Figure 6). These parameters were: the normalization constant (c_f = 0.01; Figure 6c; Table 2) of module 2 and the fraction of water routed into the fast-responding reservoir of event water (q_e = 0.63; Figure 6e), the mean transit time (MTT) of the fast fraction of event water (k_fe = 20.62 h; Figure 6f), the time delay for the event fraction response (e_lag = 0 h; Figure 6h), the fraction of water routed into the fast-responding reservoir of pre-event water (q_p = 0.27; Figure 6i), and the MTT of the fast fraction of pre-event water (k_fp = 3.84 h; Figure 6j), and the time delay for the pre-event fraction response (p_lag = 0 h; Figure 6l) of module 3. Both parameters of module 1 (S_o; Figure 6a and w; Figure 6b), the memory timescale parameter (w_f) of module 2 (Figure 6d), and the slow fractions of event (k_se; Figure 6g) and pre-event (k_sp; Figure 6k) water of module 3 showed equifinality [79]. That is, these parameters did not reach a single peak associated to the highest NSE in their parameter space distributions. The distributions of the model parameters for model structure 4 calibrated using δ¹⁸O was similar (Figure S9) and the ranges of calibrated parameters using both tracers were in strong agreement (Figure 7).

Similar results were found for model structure 2 using EC for calibration. For this structure, the module 3 parameters q_e, k_fe, e_lag, q_p, k_fp, and p_lag, as well as the constant fraction of effective rainfall routed as event water parameter f (module 2) reached a single peak in the parameter space yielding the highest NSE (i.e., seven out of 11 parameters; Figure S6). However, even tough for the same model structure calibrated using δ¹⁸O, the parameters q_p and k_fp did not reach a single peak in the parameter space (i.e., only five out of 11 parameters did; Figure S5), the ranges of calibrated parameters using both tracers were in agreement (Figure S11).

For model structure 1 calibrated using EC, none of the parameters of module 1 reached a single peak in the parameter space, whereas the parameters f (module 2), and k_e, e_lag, k_p, and p_lag (module 3) did (i.e., five out of seven parameters; Figure S4). The calibration of structure 1 using δ¹⁸O, showed that the parameter w (module 1) reached a single peak in the parameter space in addition to those that did for the model calibrated for EC (i.e., six out of seven parameters; Figure S3). Similarly, the ranges of the calibrated parameters for both tracers were not in agreement (Figure S10). For model structure 3 calibrated using EC, only the c_f parameter (module 2), and the k_e, k_p, and e_lag parameters (module 3) reached a single peak in the parameter space (i.e., four out of eight parameters; Figure S8). The same parameters reach a single peak in the parameter space when structure 3 was calibrated using δ¹⁸O (Figure S7) and the ranges of calibrated parameters using both tracers were in agreement (Figure S12).

3.5. Comparison of Model Results for Both Tracers

Even though the TTDs of the event fractions using both tracers were similar for transit times (τ) longer than 20 h (Figure 8), there was some discrepancy in the τ distributions shorter than 20 h. For shorter τ, the calibration with δ¹⁸O yielded a TTD with a predominance of shorter τ in comparison to the calibration of the model calibrated with EC (Figure 8a). The shape of the TTDs of the pre-event water fractions (Figure 8b) were opposite to those of the event water fraction. That is, for shorter MTTs (<10 h), there was a predominance of shorter MTTs for the calibration with EC compared to the calibration with δ¹⁸O.

4. Discussion

4.1. Selection of the Best Model Structure Using Water Isotopes and Electrical Conductivity for Model Calibration

The various TraSPAN model structures can be used to test different hypothesis about the processes that control streamflow generation [69,70] or as a rejectionist framework to improve knowledge about the hydrological behavior of a given catchment [26]. Thus, the modeler can hypothesize and test different hydrological behaviors and assumptions in a catchment by building and testing conceptual models that represent different catchment functions [80,81].

In our study, the four tested model structures (i.e., competing hypotheses [82,83]) provided strong fits in terms of the NSE objective function using both tracers for calibration (NSEs > 0.79; Table 3). However, further evaluation of the models’ parsimoniousness, depicted important differences among the evaluated structures (Table 3). The simplest TraSPAN structures, including a single reservoir for the event and pre-event water transit times assuming either constant (structure 1 with 7 parameters; Figure 2a) or a time variant (structure 3 with 8 parameters; Figure 2c) fraction of P_eff routed as event water, provided less parsimonious results than the more complex structures 2 and 4. Between structures 2 and 4, which represented the catchment response with two connected linear reservoirs for each the event and pre-event water fractions and a constant (structure 2 with 11 parameters; Figure 2b) or time variant (structure 4 with 12 parameters; Figure 2d) fraction of P_eff routed as event water; structure 4 was the most parsimonious to represent the internal catchment response using both tracers (Table 3).

Regarding the results of the simulations of the observed Q and tracer data, the simplest structures 1 and 3 poorly resembled the hydrograph. These structures did not reproduce the peak responses and the recessions that followed them (DOY 306–311; Figure 4a,d, Figure S2a,d). Regarding the simulation of the tracer dynamics, these structures captured relatively well the δ¹⁸O_Q decrease and EC_Q increase during the hydrograph falling limb but poorly resembled the tracer dynamics at the beginning of the event (Figure 4b,e, Figure S2b,e). These results indicate that the hydrological behavior of Mack Creek is not well represented when a single reservoir is used to route the event and pre-event water fractions regardless of the treatment of the P_eff as constant or time variant. On the other hand, structures 2 and 4 provided better representations of the hydrograph and the tracers’ dynamics. Both of these structures calibrated with both tracers were able to better resemble the highest Q peak (DOY 305) and the following rapid recession (Figure 5a,d, Figure S1a,d). These results suggest that the Mack Creek hydrological system is better represented by two connected linear reservoirs for each the event and pre-event water fractions. Regarding the tracer dynamics, structure 4 provided the best fits of the tracer data in comparison with structure 2 for both tracers (Figure 5b,e, Figure S1b,e). These results indicate that the system is not only better represented by two connected reservoirs, but that the internal mixing of the tracer is better represented when the P_eff routed as event water is treated as time-variant (structure 4) than when it is treated as constant (structure 2). Even though not directly investigated at Mack Creek evidence of non-stationary flow conditions and tracer dynamics have been shown for WS10, a 0.1 km² catchment located near the outlet of the H.J. Andrews Forest (Figure 1) [63,84,85].

4.2. Comparison of Flow Partitioning Modelling Results Calibrated for Water Isotopes and Electrical Conductivity

The ranges of the calibrated parameters for TraSPAN structure 4 were remarkably similar using both WSI and EC for calibration (Figure 6 and Figure 7, and Figure S9, and Table 2). In addition, the temporal variation (Figure 5c,f) and the estimated proportions of pre-event water yielded by the model calibrated for each tracer were in agreement (74.9% for δ¹⁸O and 81.0% for EC, Table 3). Even though previous hydrograph separation investigations using simple tracer mixing models have found contradictory results regarding the reliability of the event and pre-event water contributions when EC is used as a tracer due to its non-conservative nature [16,18,48,50,51,53,86,87,88], our findings suggest that this issue could be related to the lack of simultaneous representation of the water flow transport in addition to the tracer mixing within a given hydrological system. We were able to account for both using our TraSPAN numerical modeling approach.

Similar findings were reported by Mosquera et al. [76] in a study conducted to evaluate the baseflow transit times in a nested system of tropical alpine (páramo) catchments in South America. These authors found that the spatial variability of transit time estimates for their nested system of catchments—based on the calibration of a lumped conceptual model using δ¹⁸O—was highly correlated with the catchments’ mean yearly baseflow EC (r² = 0.89) and suggested that average EC values could be used as an inexpensive proxy to estimate baseflow transit times. Even though these results cannot be generalized until additional investigations are carried out in different environments, the results from Mosquera et al. [76] and our study highlight the potential advantages of high-resolution monitoring of EC for hydrological modelling applications, particularly when flow dynamics are accounted for.

4.3. On the Use of High-Temporal Resolution EC in Hydrograph Separation Modelling

It is worth noting that the structure 4 simulation results yielded better model performance, both in terms of goodness of fit and parsimoniousness, when calibrated with the finer temporal resolution of EC data than when calibrated for δ¹⁸O (Table 3). The simulation results also showed that a higher accuracy to represent the hydrograph and tracer dynamics was obtained when structure 4 was calibrated using the higher-resolution EC data (sub-hourly), than when calibrated for WSI data collected every 4 h. For instance, the hydrograph highest peak and following recession on DOY 305 (Figure 5a,c) and the rapid δ¹⁸O increase/EC decrease on DOY 304–305 (Figure 5b,d) were better represented by the model calibrated using EC data. These results indicate that even though the model calibrated for both tracers depicted similar hydrological behavior (see Section 4.4 for details), the calibration for the higher-resolution EC data allowed better capturing of the fast occurrence of flow transport and solute mixing processes. These simulations resulted in the release of higher amounts of pre-event water for the model calibrated for EC (81%), particularly during the highest hydrograph peak on DOY 305 (Figure 5f), with respect to the calibration for δ¹⁸O (75%; Figure 5c).

Other investigations have also reported higher accuracy to represent the temporal variability of the tracer’s signal, and further interpretation of water routing in hydrological systems when using high-temporal resolution tracer data [89,90]. Our results do not only support the findings from these authors, but also demonstrate the potential of the inexpensive collection of EC data at high-temporal resolution to improve understanding of fast occurring (seconds to minutes) hydrological processes.

4.4. Process-Based Understanding of Hydrological Behavior

We evaluated the similarity in hydrological behavior of the transport and mixing of water within the Mack Creek catchment based on the shapes of TTDs provided by the best model structure calibrated for each tracer. Despite the overall similarities of the TTDs’ shapes for both tracers (Figure 8) with dominance of long transit times (τ > 10 h), our results showed small differences for shorter transit times (τ < 20 h for the event water fractions, Figure 8a; τ < 10 h for the pre-event water fractions, Figure 8b). This discrepancy could be related to different temporal resolution of the data used for calibration, as reported in other studies [89]. In other words, the use of finer temporal resolution EC data (every 0.25 h) in comparison to the WSI data (every 4 h).

From a process-based perspective, the TTDs obtained from the model calibration using both tracers indicate that the catchment acts as a connected system of two water reservoirs each with a fast and a slow transit function contributing large amounts of pre-event water (75–81%) to discharge (Table 3). These two reservoirs likely represent: (1) fast event and pre-event water moving through the shallow permeable soils and (2) slow event and pre-event water moving through the fractured parent material (hereafter referred to as the groundwater reservoir, GW). That is, this catchment has poorly developed gravelly loam soils with high infiltration rates (>500 cm h⁻¹) that are overlaying highly weathered and fractured bedrock [57]. Baseflow MTT at Mack Creek was estimated around 2 ± 0.49 years, i.e., one of the longest within the H.J. Andrews catchments [57]. This MTT estimation indicates that water is likely stored in the GW reservoir, as has been observed in other headwater catchments with relatively permeable geology (e.g., [73,91,92]). Stored water that can be released rapidly to streams during rainfall events given the connectivity between hillslopes and riparian areas. Such connectivity has been observed near our study area at WS10 (Figure 1) [63]. Additionally, even though high amounts of pre-event water were released as streamflow during the monitored event, the total runoff (20.6 mm) accounted for only 13% of the total precipitation input (155.1 mm). Again, given the permeable soils and fractured bedrock of the catchment, and the dry antecedent moisture conditions during the monitoring period (28 October–7 November)—which corresponded to the beginning of the fall season after the summer period in the particularly dry water year 2015—only a small fraction of total precipitation was converted into P_eff. These results from the hydrometric analysis indicate that the rest of water inputs must have filled the initially low water storage of the GW reservoir after the dry summer. In this context, the strong goodness of fit of our model to simulate the high temporal variability of the streamflow and tracer data dynamics (Figure 5a–d) indicates that through the simultaneous calibration of water and solute fluxes, our model was not only capable to successfully account for the water flux and tracer mixing dynamics, but also for the recharge of water in the groundwater system.

5. Concluding Remarks and Future Directions

Our evaluation of a tracer-based hydrological model (TraSPAN) demonstrated that the assessment of different modeling structures—each representing a different rainfall-runoff response—allowed for a better identification of the hydrological system functioning. For our analyzed rainfall event at the H.J. Andrews experimental forest, we found that the same model structure was best at representing hydrographic and tracer dynamics using either electrical conductivity/specific conductance (EC) or water-stable isotopes (WSIs) collected at high temporal resolution for calibration. The model results using both tracers not only showed a remarkable agreement to fit the observed data and the calibrated parameter distributions, but also in terms of the process-based understanding of the hydrological system. Moreover, the use of sub-hourly (0.25 h) collected EC data also allowed to better simulate the catchment hydrological response in comparison to the 4 h δ¹⁸O data. These findings highlight the potential of using low-cost EC data collected at high temporal resolution in combination with a flexible hydrological modelling framework to better understand catchment hydrological behavior.

Despite the advantages of our applied methodology, we acknowledge that before the widespread application of our approach, future research should focus on understanding site-dependent geochemical conditions under which the applicability of EC as a conservative tracer is suitable. Future research should also focus on implementing monitoring strategies that allow for the combined collection of water isotopes and geochemical data at the highest possible temporal resolution during a variety of climatic conditions (i.e., wet, transition, and dry periods). We encourage such efforts, as these data would provide information about the spatial and temporal variability in weathering rates to better establish the conditions under which the use of high-resolution EC observations would yield robust datasets to investigate important, but yet unresolved, questions in catchment hydrology that require the identification of fast-occurring hydrological processes, while reducing monitoring costs and the degree of uncertainty in simulated flow components.