Stream ﬂ ow and hydrogen ion interrelationships identi ﬁ ed using data-based mechanistic modelling of high frequency observations through contiguous storms

With the aim of quantifying the purely hydrological control on fast water quality dynamics, a modelling approach was used to identify the structure (and dynamic response characteristics or DRCs) of the relationship between rainfall and hydrogen ion (H þ ) load, with reference to rainfall to stream ﬂ ow response. Unlike most hydrochemistry studies, the method used makes no a priori assumptions about the complexity of the dynamics (e.g., number of ﬂ ow-paths), but instead uses objective statistical methods to de ﬁ ne these (together with uncertainty analysis). The robust models identi ﬁ ed are based on continuous-time transfer functions and demonstrate high simulation ef ﬁ ciency with a constrained uncertaintyallowinghydrologicalinterpretation ofdominant ﬂ ow-paths and behaviour of H þ loadin four upland headwaters. Identi ﬁ ed models demonstrated that the short-term dynamics in H þ concentration were closely associated with the stream ﬂ ow response, suggesting a dominant hydrological control. The second-order structure identi ﬁ ed for the rainfall to stream ﬂ ow response was also seen as the optimal modelforrainfalltoH þ load,evengiventheverydynamicconcentrationresponse,possiblyindicatingthe same two ﬂ ow-paths being responsible for both integrated responses. Within this study the latest routines for identi ﬁ cation to simi-larly cant advances in the for continuous-time (Young & ; & is with the parallel in


INTRODUCTION
Short-term pulses of stream acidity can have a disproportionate negative effect on stream ecology (Lepori et al. ; Monteith et al. ), but models specifically of these short-term acid pulses in storms have received little attention since the early 1990s. Further, it is increasingly recognised that water quality sampling and/or in situ monitoring should be undertaken continuously at a high enough frequency to fully capture the dynamic behaviour of solute concentrations and load within streams (Kirchner & Neal ). Studies using continuous higher frequency sampling or in situ monitoring of a wider range of stream water quality variables are increasing (Table 1). However, even some of these studies may have insufficient sampling to capture the rapidly changing concentrations through storms within small headwater streams.
Several dynamic models capable of simulating changing stream acidity through storms were developed in the 1980s.
These included the Birkenes model (Christophersen et al. ), MAGIC model (Cosby et al. ) and CAPTAIN model (Langan & Whitehead ). Significant differences in those catchment processes that are assumed to be the dominant controls of stream acidity exist between these models (Christophersen et al. ). The original MAGIC model of Cosby et al. () assumed that the time series of stream H þ concentration (with a weekly simulation time-step) is primarily a function of a time-varying geochemistry of a single soil unit combined with the effects of CO 2 de-gasing as soil-water exfiltrates into a stream. In contrast, the Birkenes model of Christophersen et al. () assumed that stream acidity is also a function of the mixing of This is an Open Access article distributed under the terms of the Creative chemically distinct soil-waters sourced from two separate soil horizons, and the implicit effects of hydrological routing to and from these horizons. Other models assuming a key role of hydrological routing or pathways do, however differ in that they assume that fewer (i.e., one: Whitehead et al. ) or greater (i.e., three or more: Grip et al. ) numbers of water pathways are required to simulate the dynamics in stream acidity. Even after 20 years of further research into water pathways controlling stream chemistry, there is often no consensus as to the geochemical and hydrological structure of catchment models (e.g., whether they include one or more dominant strata-specific pathways or separate macropore flows) needed to simulate stream acidity or other water quality variables (Kirchner ). Advances in our understanding of uncertainty in hydrochemical modelling (e.g., Medici et al. ) do, however, provide ever greater evidence that with overly complex models many different combinations of processes (hydrological, geochemical or both) can simulate stream acidity. Consequently, accurate simulation of an output variable such as stream H þ concentration does not necessarily mean that the dominant mechanism for the release of H þ to streams has been represented accurately. Correct outputs may be simulated for the wrong reasons (Beven & Westerburg ). Models that have many more parameters than can be justified by the complexity of the dynamics of the variable to be simulated (i.e., overly complex models) are now seen as poor models because of the reduced identifiability (i.e., increased parameter uncertainty) arising from the routines used to identify optimal parameter sets. As a result of these two issues there is an increasing acknowledgement of the value of using parsimonious (i.e., efficient but simple) model structures (i.e., those that do not include too many hydrological or geochemical processes or pathways). Similarly, there is an awareness of the value of not fixing a model structure without first evaluating alternatives (e.g., not forcing the simulation of acid runoff to follow two water paths without also assessing the simulation capabilities of a single or numerous water pathways). This issue was appreciated following the early hydrochemical modelling work (e.g., Christophersen et al. ) but often ignored subsequently.
The philosophy of the data-based mechanistic (DBM) modelling approach adheres to both of these principles (Young ). The first stage in the DBM approach is the application of a large number of mathematical relationships (often in the form of transfer functions), to describe the dynamics between observed input time series (e.g., rainfall, temperature), and output time series (e.g., H þ concentration), and thereby form data-based models. Those identified models that do not meet strict mathematical-statistical
Between 1957 and 1959, LI3 was almost completely afforested with a combination of Sitka spruce (occupying 53% of maximum forest cover), Lodgepole pine (16%), Japanese larch (13%) and a Lodgepole pine/Japanese larch mix   Measurement of pH was undertaken using digital differential pH sensors (model: pHD-SC) combined with an SC200 controller (Hach Lange, Düsseldorf) and recorded at 15-minute intervals using further CR1000 data-loggers.
The sensor was recalibrated at intervals determined by the sensor-controller system. In situ monitoring of other water quality variables (turbidity, nitrate, DOC, total organic carbon, colour and EC) was also undertaken at the locations of the pH sensors, but these data are not reported here.
All data were downloaded from the data-loggers on a weekly basis, and manual pH and water depth readings used to check for sensor malfunction; no such malfunctions were observed. Data were quality assured using dedicated CT transfer function modelling using the DBM approach The first stage of the DBM modelling process involved fitting a range of transfer function models to the data. For production of model parameter estimates in highly dynamic input, single output system can be given as follows: where u(t) is the deterministic input signal, x(t) is the noisefree output signal, y(t) is the noisy output signal, ξ(t) is the noise and τ is the pure time delay (i.e., delay between an input and the first response in an output) in units of time.
A(s) and B(s) are polynomials in the derivative operator s ¼ d/dt as follows: where n and m can be any positive integers with m n.
These equations are able to be expanded for multivariable systems, i.e., multiple input (e.g., rainfall and temperature) and single output (e.g., streamflow) or multiple output (e.g., streamflow and H þ ).
When input, u(t), and output, x(t), observations have been made at a constant sampling interval, Ts (i.e., in discrete time), then the sampled signals are represented as u(t k ) and x(t k ) respectively, and the output observation equation is as follows: where x(t k ) is the sampled value of the noise-free unobserved output (x t ), which for data sampled in DT is assumed to be contaminated by noise, ξ(t). The noise is inde- Within this routine, high-order noise is removed from the data (i.e., dynamics in the system that are much quicker than the quickest time constant) and which hinder identification of true model structure. The skill in the RIVCBJID routine also arises from its use of information within the covariance matrix (see Box et al. ).
Model structure is typically presented as a triad of denomi- The identified models are assessed statistically using a measure of simulation efficiency R t 2 and the Young Information Criterion (YIC). The R t 2 is given as Equation (5) where σ 2 error is the variance in the model residuals and σ 2 obs the variance in the observed data. This measure of simulation efficiency is equivalent to the simplified form of the Nash and Sutcliffe statistic and so differs from the R 2 that uses the variance in the model output as the denominator.
The YIC is an objective statistical measure given as: The first term within Equation (6) is a measure of the model efficiency (see Equation (5)) and the second term, normalised error variance norm, is a measure of the degree of over-parameterisation (Young ). Generally, as model complexity increases, so does the ability to capture more and more of the dynamics in the observed time series, but at the expense of increasing parametric uncertainty.
There is, therefore, an optimum complexity (or 'model order') given that more complexity increases the uncertainty in the individual parameters identified. The YIC is a measure of whether the model has become overly complex (i.e., too many parameters) given the amount of information contained within the observed data series. A change of þ1.0 or greater in YIC as model order is increased indicates that the model structure has become too complex for the information contained in the time series, i.e., overparameterised (Ockenden & Chappell ), and the simpler model should be accepted as the optimal order.
The procedure for identifying the optimal model structure using these two measures begins by identifying the first-order model (from those convergent models) with the highest R t 2 that does not have complex roots (see Box Where second-order models have a higher R t 2 than that of the optimal first-order model, these models are then examined for a change in the magnitude of the YIC from first to second-order models. If a second-order model: (i) has a higher R t 2 than the optimal first-order model, (ii) its YIC is less than þ1.0 different to the firstorder model, (iii) it does not have complex roots, (iv) it does not exhibit oscillatory behaviour, and (v) the sign (±) of the two identified roots are the same, then this model is now accepted as the optimal structure. This procedure is repeated for third-order models, then fourth and so on.
Models with a structure higher than sixth order are not normally examined because they are invariably overparameterised. Typically, with rainfall to streamflow models attempts are made to identify models with between zero and six pure time delays. Consequently, by attempting to identify models from first-order to sixth-order with zero to six pure time delays, this gives 182 possible model structures to be identified using RIVCBJID. Many of these possible structures will not give convergent solutions.
Within this study, uncertainty in the model parameters was assessed using 1,000 Monte Carlo realisations of the optimal A and B terms (Equations (2) and (3) where ▵t is time-step of the observations (i.e., 15 min), SSG 1 is the SSG of the slow response pathway and SSG 2 is the SSG of the fast response pathway. The TC is the residence time of the response (i.e., not the residence time of a water particle or H þ atom) between a rainfall input and either a A range of purely linear models from first-order structures to sixth-order structures were identified, and each had a pure time delay between the onset of a rainfall event and a response in the stream of between zero and 90 minutes (i.e., six observation time-steps). The most efficient 20 models identified (ordered by R t 2 ) for one of the four catchments (LI7) as an example, are shown in Table 4. The optimal models (i.e., high R t 2 , but without a ! þ1.0 change in YIC compared to lower order models) for all four catchments are shown in Table 5.
The use of purely linear CT-TFs has produced optimal models with at least 90% efficiency (R t 2 ) of streamflow simulation (   The term den is the number of transfer function denominators (recession or 'a' parameters), num are the number of transfer function numerators (gain or 'b' parameters), delay is the pure time delay between rainfall and runoff response, YIC is the Young Information Criterion and Rt 2 is the efficiency measure. Terms BIC (Bayesian Information Criterion), S2 and condP are additional measures of model over-parameterisation used by other studies, so included here for reference. The model considered to be optimal (i.e., high R t 2 , but without an integer loss of YIC compared to lower order models) is highlighted in bold and underlined. Model structure is given as [Denominators, Numerators, Pure Time Delays]. The fast% term is percentage of response following a fast pathway of a second-order model (see Equation (9)), while the term SD of fast TC is the standard deviation of the fast TC from 1,000 Monte Carlo realisations. The steady state gain (SSG) is not presented as the rainfall time series used is an arithmetic average of the 15-minute totals for LI3 and LI6 rain gauges, rather than catchment-average totals per 15-minute intervals. whereas other forms of decomposition (e.g., feedback systems) lack a feasible hydrological interpretation. This parallel model comprises two TCs and the proportion of response associated with each TC (Table 5). The uncertainties in the TCs (expressed as standard deviations derived by Monte Carlo analysis) are also presented in Table 5.
For the LI3 results, the frequency distributions of the fast and slow TCs are shown in Figure 3 as an illustration. A narrow range of uncertainty for both the slow and fast TCs was observed from this Monte Carlo analysis, with the exception of the slow response component for LI6 (Table 5).
A mechanistic hydrological interpretation of the identified TCs can be achieved with reference to such values derived by similar modelling approaches previously applied to a range of sites each with different streamflow generation mechanisms ( Table 6). The time series of the simulated fast component for each of the four basins is shown in Figure 4(a) and the slow component in Figure 4(b).
The TCs for the fast component for the four experimental basins at Llyn Brianne ranged from 2.40 to 5.53 hours (Table 5). These values are, for example, much longer than those of measured overland flow on a tropical hillslope at 5 min or 0.08 hours (Chappell et al. ), but comparable to the response of a hillslope system elsewhere in the Cambrian mountains (i.e., 2.9 hours:   To determine the transferability of the second-order linear models of rainfall-streamflow to a later winter period that also has H þ data, the RIVCBJID algorithm was applied independently to a series of three storms that occurred between 5 and 18 February 2013. Again, a range of models with firstto fourth-order structures were identified with a pure time delay between zero and six time-steps. The optimal models and associated parameters for the four catchments are shown in Table 5 alongside the results for the earlier period. For this February period, purely linear second-order models [2 2 τ] were again identified as the optimal models and produced R t 2 efficiency values in excess of 85% (Table 5; Figure 5). This provides conditional validation (Young ) that rainfall-streamflow response at Llyn Brianne during the winter consistently has a second-order model structure. For this second period the TC for the slow component doubled, however, this was offset by a significant reduction in the proportion of the response along the slow pathway ( Table 5). The differences in the pure time delays between the two periods are very small, changing by zero to only two time-steps. Furthermore, there is no systematic shift in the fast TC between the two periods (Table 5).

Streamflow to H þ concentration modelling
Given Purely linear first-order CT-TF models, the simplest form of a dynamic model that may be presented, were able to identify between 82 and 97% of the variance in the H þ concentration for the four streams (Table 7, Figure 6). For linear second-order models, the explanation ranged from 93 to 99% (Table 8; Figure 7). The improvement in model efficiency was most marked with LI3 where it increased from 82 to 93%. Furthermore, for LI3 the second-order structure better captured the later recessions (Figure 7), and had a   Table 7). The variability between the four streams in their H þ concentration response to streamflow (Table 7) is considerably larger than that between rainfall and the fast component of streamflow (Table 5); the coefficient of variation in the TCs of the former is 18% and the latter 91% (derived from data in Tables 5 and 7). This greater variability   (Table 5) which has over 95% forest cover of about 40-year-old conifers.

Streamflow to H þ load modelling
The DBM modelling of the rainfall to streamflow relationship indicated that the response could be described by two components with quite distinct responses (i.e., TCs). Given that the H þ concentration appears to be closely associated with the hydrological dynamics in the streamflow time series at Llyn Brianne (e.g., Figure 6), the investigation next addressed whether the export of H þ (i.e., load) was associated with the two distinct response pathways identified for the streamflow generation. The final stage in this modelling was the identification of the relationship between rainfall and H þ load, and how many separate components the time series of H þ load may be divided into. If the load model has the same model structure as inferred for the streamflow response to rainfall then this allows comparison of the characteristics of the two components of the response with those of the streamflow generation pathways (i.e., the TC and the portion of response along each pathway).   As with the rainfall-streamflow modelling a range of model structures from firstto sixth-order were investigated. As with the hydrological model, purely linear second-order models were optimal for all four sites (Table 9; Figure 10). This time series of H þ load could then be readily decomposed into  The second-order models do explain two-thirds of the variance in the H þ export (i.e., R t 2 0.711 to 0.752) and are, therefore, models amenable to physical interpretation. The reduction in efficiency of modelling H þ export rather than streamflow from rainfall was to be expected given the noise observable within the H þ records, and the storm to storm variability in the H þ behaviour observed previously at Llyn Brianne ( is also observed within raw observations, and underlines the importance of high-frequency sampling to capture these differences between the H þ response on the rising versus falling stages of stream hydrographs. Where this hysteresis is present, this means that simple relationships with no dynamic storage effects between streamflow and H þ export would be highly uncertain and should be avoided for estimation of long-term loads (Littlewood ).
The TC of the fast component of H þ export for LI3, LI6 and LI7 is similar at 2.5, 2.3 and 3.4 hours, respectively; however, the export from LI8 has a TC of 7.5 hours (Table 9). As with the observation of the H þ concentration     Total masses of H þ exported were very different between the grassland and forest basins (Table 9); however, the TCs for the H þ load were very similar for both the grassland and forest basins, indicating the importance of the basin hydrology over the impact of vegetation on biogeochemical dynamics.

Minimum monitoring interval
The Nyquist-Shannon sampling theorem indicates that monitoring or sampling must be at least twice that of the fastest TC, and in practice, at least six times faster (Young ). Table 10 shows the minimum sampling interval based upon all fast TCs for the rainfall to H þ load (from Table 9), and for reference the rainfall to streamflow (from  (Table 10) The chemical concentration or load response of a stream during a storm (e.g., Figure 7) can be called the storm chemograph. An estimate of the TC of the recession of such a storm chemograph can be gained by finding the time from the chemograph peak to 63% of the return to the pre-storm concentration or load, assuming an exponential recession from a linear store. Table 11 shows the

CONCLUSIONS
Eight novel conclusions have arisen from this work: 1. The hydrological control of dynamics in the constituent H þ concentration was affirmed from the strength of the dynamic relationship with the hydrological measure of streamflow (R t 2 > 89%), i.e., the lumped measure of often complex soil-water pathways, and potentially H þ pathways, within headwater catchments.
2. The H þ load modelling indicated that two separate response components were present within the relationship (i.e., second-order dynamics). This clear bimodal behaviour was also present within the optimal rainfall to streamflow models, suggesting that H þ load was strongly moderated by the effects of these two hydrological response pathways. These two response components were identified from the observations by application of two objective statistical criteria (R t 2 , YIC), rather than assuming a priori a certain model structure, as in most non-DBM modelling studies.
3. The minimum sampling interval identified by the transfer function modelling indicates that for simulation of H þ load requires minimum sampling intervals of 23, 25, 34 and 75 min for the headwater basins LI6, LI3, LI7 and LI8, respectively. These rates are within the range estimated from a preliminary set of key studies with high frequency sampling of chemical concentration through storms. Most water quality data-sets, however, lack such data (even for individual storms), and so may be inappropriate for developing or evaluating models that can capture the hydrological controls on water quality dynamics.
4. Comparison of the identified TCs of the hydrological response characteristics with those identified within other basins having independent corroboration of water pathways, indicated that a soil-water pathway may be responsible for the fast pathways at the Llyn Brianne basins, while the slower path (during the wet winter conditions) may be associated with pathways in underlying drift strata. While other studies and many temperate sites have suggested the importance of these two pathways, they have not been previously identified by inter-comparison of DRCs derived for reference sites (Table 6).
Independent evidence of these two pathways within the studied Brianne basins is, however, required to strengthen or correct these tentative conclusions. A new measurement programme has been initiated to obtain direct observations of the two pathways (rather than estimate from their resultant impact on the streamflow response).
5. The response characteristics of the slower pathway are the least certain given the short duration of the modelled record relative to the TC of the slower pathway. Simulation of longer records is required to constrain the uncertainty within this TC. Furthermore, longer-term simulation that also includes drier conditions may allow the identification of even slower components; perhaps via hydrologically active rock fractures known to play a role within the hydrochemistry elsewhere within the Lower Palaeozoic Cambrian Mountains.
6. Relative to the purely hydrological response, the response of the H þ load was considerably faster, indicating very rapid mixing and/or ion exchange processes to an input of rainwater to these catchments. Both fast and slow components of the load response did, however, show clockwise hysteresis relative to the hydrological response components indicating a supply-limited process or exhaustion of the H þ ions along the inferred pathways with the progression of a storm. Again, these pathway-specific interpretations need to be evaluated further following the planned new programme of continuous monitoring of H þ concentration at different depths (i.e., overland flow, soilwater, water in the glacial till and water in rock fractures).
7. The basin that is primarily covered by mature conifers