BaHys—A Bayesian Modeling Framework for Long‐Term Concentration‐Discharge Hysteresis: A Case Study on Chloride

Improving river water quality requires a thorough understanding of the relationship between constituent concentration and water discharge during runoff events (i.e., C‐Q hysteresis), which may be strongly non‐linear. Analysis of C‐Q hysteresis on large temporal scales provides unprecedented insights into event dynamics and long‐term concentration trends in surface and groundwater. Despite the increasing availability of time series data on water quality, there are still limited quantitative modeling frameworks that enable this analysis. Here, we combine Bayesian modeling and an existing mass balance to model long‐term C‐Q hysteresis dynamics in multi‐decade constituent concentration and water discharge time series. We focus on the case study of chloride and demonstrate that our model can simultaneously characterize the size and rotation of C‐Q hysteresis, and diffuse and low‐flow inputs to constituent concentration using only time‐series data from the river Rhine. Over 28 years of monitoring, we find that chloride exhibits a dominant clockwise dilution behavior that does not vary considerably under different hydro‐climatic conditions, hinting to similar mobilization mechanisms over time. We also show decreasing chloride concentrations in surface and groundwater, due to the cessation of mining activities in the Rhine. Our approach uses uncertainty estimates to show the range within which model parameter values lie, aiding decision‐makers in a robust assessment of river water quality. We conclude that Bayesian modeling of C‐Q hysteresis provides a powerful framework for investigating long‐term contamination dynamics that can be extended to several constituents to find factors controlling their export, ultimately suggesting mitigation measures for river contamination.


Introduction
River ecosystems are under continuous threat from human activities (Borgwardt et al., 2019;Palmer et al., 2009; C. J. Williams et al., 2016).Activities like intensive agriculture release a variety of constituents-solutes and nutrients such as chloride, nitrogen and phosphorus-in rivers, increasing their concentration in surface and groundwater and endangering life below and above water (Garnier et al., 2016;Kaandorp et al., 2021;Oenema et al., 2005;Winter et al., 2021).Increased salinity due to chloride concentrations above the permissible limits changes the physicochemical properties of river water, threatens aquatic biodiversity, and corrodes drinking water pipelines (Dugan & Arnott, 2022;Pieper et al., 2018;Szklarek et al., 2022).High nitrogen and phosphorus concentrations degrade water quality through, for example, acidification and eutrophication (Conley et al., 2009;Krusche et al., 2003;Wang et al., 2020).Reducing such pressure requires understanding the dynamics and behavior of these constituents in rivers, including mixing mechanisms, trends, and concentration-discharge timelags (Kaandorp et al., 2021).This enables water authorities to adhere to the goals set forth by the Sustainable Development Goal 6 and by the European Water Framework Directive (EU WFD), landmark among water policies demanding control and reduction of pollution toward good freshwater status and sustainable use (Brack et al., 2017;Krueger et al., 2009).
The relationship between constituent concentration (C) and water discharge (Q), commonly known as C-Q relationship, is valuable for understanding the behavior of critical constituents in rivers (Godsey et al., 2019).This enables the quantification of reactive transport and mixing phenomena, emission, mobilization and export behavior (Guo et al., 2022;Knapp et al., 2020;Minaudo et al., 2019;Wymore et al., 2019).Analyzing the C-Q relationship across different time scales offers several insights.Long-term analyses are crucial for modeling constituent transport processes and characterizing water quality trends and seasonal dynamics in a catchment (Hirsch et al., 2010;Winter et al., 2021).However, analyses on the event scale are particularly relevant to describe the non-linear C-Q relationship, which can reveal mixing and emission behaviors of the constituent of interest (Campbell et al., 2005;House & Warwick, 1998;Lloyd et al., 2016b).This non-linearity results in the phenomenon of hysteresis.
Hysteresis indicates a dependence of constituent concentration on previous concentration values, as well as on instantaneous water discharge (Gonzalez-Nicolas et al., 2021).This phenomenon frequently occurs during runoff events, where an increase in water discharge caused by short-lived hydrologic events such as storms, snow melts, or irrigation, corresponds to a higher constituent mobilization (Burns et al., 2019).Hysteresis originates due to multi-factorial hydrochemical mechanisms, including delayed or anticipated emission of constituents into the stream or low-and high-frequency sources of runoff mixing at different times (Bieroza & Heathwaite, 2015;Musolff et al., 2021).These mechanisms ultimately reflect the impact of internal and external stressors, such as land use, human activity, hydrochemical and climatic conditions, on constituent concentration (Burns et al., 2019).
Hysteresis is commonly evaluated through hysteretic loops, resulting from the direct comparison of each runoff event and constituent concentration.The loop magnitude, shape and rotation characterize the constituent export regime and the proximity of its source to the measurement point (Rose et al., 2018;Zuecco et al., 2016).Such an analysis is often performed manually, but it becomes time-consuming and burdensome for the increasingly larger quantities of water quality data that are becoming available, requiring automated evaluation for robust and comparative interpretation (Millar et al., 2022).To support such automation, various hysteresis indices have been developed over the past decade (Butturini et al., 2008;Lawler et al., 2006;Lloyd et al., 2016a;Vaughan et al., 2017;Zuecco et al., 2016).For example, Lloyd et al. (2016a) quantified the size and rotation of hysteretic loops by using differences in concentration between the rising and falling limbs of the normalized runoff event (Lloyd et al., 2016a), Zuecco et al. (2016) relied on definite integrals for this estimation.To study the C-Q relationship, numerical models have also been developed.These models are often based on a power law, C = aQ b , where the exponent b diagnoses the chemostatic or chemodynamic behavior of a catchment (Basu et al., 2010;Godsey et al., 2009;Jawitz & Mitchell, 2011;Musolff et al., 2015;Seibert et al., 2009;Thompson et al., 2011).While this law has been used effectively to describe the C-Q relationship at different temporal scales (Basu et al., 2010;Minaudo et al., 2019), it fails to model hysteresis (Musolff et al., 2021).Models integrating hysteresis into the power law for individual runoff events have proven effective in describing the mechanisms governing hysteresis in rivers (Bowes et al., 2005;House & Warwick, 1998;Minaudo et al., 2017), including the mixing of diffuse and low-flow inputs (Bowes et al., 2005;House & Warwick, 1998).However, the extension of these models over the long term has only recently been researched.For example, Dupas et al. (2018) studied longterm hysteresis to explore multi-decade nitrogen and phosphorus dynamics; Musolff et al. (2021) recently showed that integrating hysteresis into the power law allows simultaneous estimation of inter-annual and event-scale phenomena for, for example, specific conductance and nitrate.The integration of hysteresis in the long term with numerical modeling allows knowledge from long-term and event-scale analyses to be simultaneously combined in a unified framework.
Here, we propose a model that combines direct observations and a mass balance to study C-Q hysteresis in the long term, including event dynamics and long-term trends of constituent concentration in surface and groundwater.Our proposed approach extends an already existing event-scale model (House & Warwick, 1998) to the constituent time series to simultaneously characterize the size and rotation of hysteretic loops, and diffuse and low-flow inputs in the long term.This study aims to assess that (a) hysteresis, including its size, rotation and slope, can be reliably reproduced for all the runoff events under temporal and hydrological variability, and that (b) diffuse and low-flow inputs can be derived from available observations alone while retaining their hydrochemical significance.To this end, we implement a Bayesian model that extends the semi-empirical mass balance developed by House and Warwick (1998) to multi-decade time series of constituent concentration and water discharge.Modeling long-term hysteresis is often limited by data availability.For chloride, a wealth of daily data is available for the river Rhine, as its monitoring is essential to ensure good water quality (IKSR, 2021).Furthermore, while hysteresis has been widely studied for nutrients such as nitrate and total phosphorus (Bieroza & Heathwaite, 2015;Bende-Michl et al., 2013;Outram et al., 2014), there are limited studies on chloride, especially in large river systems such as the Rhine.Therefore, chloride is a relevant case study to show the ability and potential of our approach to model long-term hysteresis dynamics.We utilize probability estimates, existing knowledge and observations to statistically assess the model's ability to describe hysteresis.We discuss the insights the model currently provides and those that could be obtained through further analyses.Finally, we examine the opportunities and challenges of this approach to other data sets with different sampling frequencies and constituents.

A Bayesian Approach to Model C-Q Hysteresis Over Time
Figure 1 presents a schematic illustration of BaHys, the approach we propose, described in detail in the remainder of this section.BaHys builds on the semi-empirical mass balance equation proposed by House and Warwick (1998) to model a single hysteretic event, presented in Section 2.1.1 (step 1).This formulation allows C-Q hysteretic behavior to be interpreted with relevant hydrochemical properties; in the second step (Section 2.1.2),we reformulate the mass balance on continuous time series to model their long-term dynamics.We assume that the runoff events are delineated in advance; Section 2.4.1 provides further details on how this was done in our analysis.This new formulation describes hysteretic dynamics with four time-varying unknown hydrochemical parameters, corresponding to constituent load at low flow (l), constituent diffuse concentration (h), water discharge at low flow (q l ) and size and rotation of hysteretic loops (p).We employ Bayesian modeling to estimate all the parameters with relative uncertainties only from observations (step 3).After defining prior probabilities for the parameters (Section 2.1.4),Bayesian modeling allows for their estimation by sampling from their posterior distributions (Section 2.1.5).After assessing model performances (step 4, Section 2.2), the entire problem formulation enables exploring the hydrochemical mechanisms governing C-Q hysteresis in the long term (step 5, Section 2.3) for the constituent of interest.The employed mass balance can be generalized to several constituents: we first provide a general description of our approach, and then show its performance on daily chloride.In the discussion, we evaluate how the framework performs on a higher-frequency data set.

Mass Balance Definition for a Single Hysteretic Event
To model C-Q hysteresis in a single runoff event, we follow the semi-empirical mass balance equation formulated by House and Warwick (1998) as: where C and Q are respectively the observed constituent concentration and water discharge defined over time t, in a single event.The equation parameters are hydrochemical properties of hysteresis; we rely on the explanation provided by House and Warwick (1998) for their definition.The constituent load at any point is the result of four main contributions: (a) inputs from point sources, such as effluents from sewage treatment plants; (b) inputs from diffuse sources, such as agricultural runoff; (c) base flow components, comprising groundwater and delayed subsurface runoff; and (d) internal processes that may increase or reduce the constituent mass due to, for example, interactions with suspended material (House & Warwick, 1998).In low-flow conditions, diffuse inputs are small, and the water composition is primarily determined by point sources and groundwater inputs.These contributions are difficult to separate and can be grouped into constituent and discharge contributions at low flow, c l and q l in Water Resources Research 10.1029/2023WR035427 CAIROLI ET AL.
Equation 1, when the diffuse contributions are minimum (House & Warwick, 1998).m net comprises the constituent mass gain or loss caused by internal processes; its contribution is difficult to measure and interpret without additional specifications (House & Warwick, 1998): to simplify the mass balance and focus on the interpretable hysteretic contributions, we consider m net as part of the constituent load at low flow, introducing the parameter l: l = q l c l + m net . (2) Diffuse contributions become particularly relevant during runoff events, hence for hysteresis, and, depending on the nature of the constituent, can result in increased constituent concentration, or constituent dilution in base flow  1, presented in the Methods section, and describes the relationship between constituent concentration and water discharge in a single hysteretic event j, depending on the parameters θ = {p, h, l, q l }.The equation in step 2 is a simplified form of Equation 4, which defines the relationship between constituent concentration and water discharge for the entire time series, according to the set of parameters Θ = {θ 1 , …, θ J }, which vary between hysteretic events.A more detailed explanation of parameters and equations is reported in the Methods section of the manuscript.
Water Resources Research 10.1029/2023WR035427 CAIROLI ET AL. (House & Warwick, 1998).In Equation 1, these contributions are described by the parameters h and p. Hysteresis is implemented through the dependence on the discharge rate dQ(t) dt : the parameter p regulates the rotation and the size of the hysteretic loop, while h characterizes the constituent diffuse concentration (House & Warwick, 1998).A direct comparison of h and the flow-weighted constituent concentration C f allows us to understand better the contribution of diffuse sources to the total constituent concentration in the river (House & Warwick, 1998).For a single event: We refer to the described hydrochemical properties as θ = {p, h, l, q l }; to model these properties on a large temporal scale, our approach adds two main contributions to the approach proposed by House and Warwick (1998): (a) we embed the equation developed for a single hysteretic event in a numerical framework that allows its extension to long continuous time series (Section 2.1.2);(b) instead of utilizing a least-squares approach, which does not allow estimation of all unknown parameters (House & Warwick, 1998), we utilize Bayesian modeling to do this directly from the observations (Section 2.1.3).

Extension of Mass Balance to Continuous Time Series
Let the constituent concentration and water discharge be continuous time series defined by N time points, and J the total number of runoff events.To extend the application of Equation 1 to J hysteretic events on the time series, we re-write the equation as: where p ∈ R N×1 , h ∈ R N×1 , l ∈ R N×1 , q l ∈ R N×1 are vectors of length N, holding time-varying parameter values (p(τ), h(τ), l(τ), q l (τ)), discretized at some time discretization step τ ∈ R, which in our experiment is one day.
Although the hysteretic patterns vary throughout the events, the hysteretic contributions may be assumed constant in each event (Musolff et al., 2021).To ensure that the parameter vectors retain such properties, we perform the following mathematical operations.Let each runoff event have a starting point s j and an ending point s j+1 .For each event, define a vector of changing points a j = [a j (1),…,a j (N)] T ∈ R 1×N , such that: Define the matrix A = [a 1 ,…,a J ] ∈ R J×N , holding the changing points for all the events, and define, for each parameter, the vector holding the parameter values for all events: According to this, define the time-varying vectors holding the parameter values in each time point τ of the time series: Water Resources Research In Equations 10-13, we define each parameter vector with constant values within each runoff event; such values change as soon as a new event starts.The new set of parameters, comprising of the parameter values for each event, can be denoted as Θ = {θ 1 , …, θ J }.Note that while the assumption of constant values within events holds for p j , which summarizes relevant features of hysteretic loops (Bieroza & Heathwaite, 2015;Musolff et al., 2021), it may be too approximate for diffuse and low-flow contributions, which can vary considerably within each event depending on its duration and the data frequency.The low-flow and diffuse estimates should therefore generally be regarded as approximate measures that capture the average diffuse and low-flow inputs in each event and can be used to evaluate long-term temporal trends and variations.BaHys assumes that the starting and ending points of each runoff event are selected beforehand, which can be done manually or with an automated algorithm (Millar et al., 2022).Both implementations are possible with our approach; which to use, depends on the goal of the study and the knowledge about the system.

Parameters Estimation
Bayesian modeling allows us to estimate all the unknown parameters (θ j ) from the observations.According to the Bayes theorem, the model learns from the parameters' prior probabilities P(θ j ) and the likelihood to assign a posterior probability P(θ j |D) to each parameter, given the observations D. Sections 2.1.4and 2.1.5report details on the definition of priors, likelihood and posterior distributions.The reader is referred to textbooks for further details on principles of Bayesian modeling (Congdon, 2007;Gelman et al., 1995;McElreath, 2020).

Definition of Priors and Likelihood
According to Equations 6-9, we define the priors for each parameter in each j event.The priors summarize the available knowledge on the unknown parameters, and should be properly selected to sufficiently constrain the parameters to plausible values (Gelman & Hennig, 2017).To this end, we utilize available contextual knowledge to specify the priors, as follows: The choice of the distribution and related parameters in the priors is specific to the studied problem; here we clarify the rationale behind these choices to facilitate the application of the proposed framework to various data sets and conditions.
1. p j : determines the rotation of each hysteretic loop and influences the size.Positive values of p j define clockwise loops, and negative values anti-clockwise loops.In the absence of hysteresis, p j is close to 0. On large time series, knowing a priori which runoff event would generate a hysteretic loop is not always feasible, especially in the absence of a priori information.However, if we assume that p j follows an appropriate probability distribution, we will only allow p j to deviate from zero if there is hysteresis.The Laplacian distribution has a faster tendency to zero rather than other continuous distributions, such as the Gaussian: in this way, we make p j diverge from zero only if the event defines a hysteretic loop.We consider a skewed distribution because hysteretic events are often asymmetric (Bieroza & Heathwaite, 2015;Liu et al., 2022; M. R. Williams et al., 2018).The distribution is defined by location μ, centered at 0 because p j should diverge from 0 only if hysteresis is present; scale σ, set at 5 because the size of the loops is not known in advance, and in this way we let the model exploring a sufficiently large probability space to find the parameter values; and skew parameter τ, set to 0.35 (i.e., right-skewed) because we expect a prevalence of clockwise loops from preliminary observations of the data.

Water Resources Research
10.1029/2023WR035427 2. h j : determines the extent of diffuse inputs.h j can only assume positive continuous values, and therefore a gamma distribution was chosen.The distribution is defined by shape α, set to 2 and an inverse scale β, set to 1.These settings allow us to obtain a finite distribution with a peak much above zero and were chosen because we expect high values for h j , related to a prevalence of diffuse versus low-flow contributions during hysteresis.3. l j , q lj : determine low-flow contributions of constituent load and water discharge.Likewise h j , these parameters can only assume positive continuous values, and therefore a gamma distribution was chosen.We define the distribution by shape α, set to 1.5, and inverse scale β, set to 2. These settings define a finite distribution with a peak above zero, but lower than h j , and a lower degree of variation compared to h j .In fact, we expect lower and more constant values for low-flow contributions compared to diffuse inputs between runoff events.
These choices resulted in weakly informative priors for our problem, meaning that the chosen distributions are spread out enough to explore a sufficiently large space to find appropriate parameter values according to the observations and the likelihood.The likelihood is the probability distribution associated with the observations, which are treated as fixed contributions, while the unknown parameters vary according to the priors (van de Schoot et al., 2021).We assume the likelihood to follow a Gaussian distribution of scale σ, after scaling the observations by their own standard deviation to ensure non-negativity:

Estimation of Posterior Distributions
The utilization of the Hamiltonian Monte Carlo (HMC (Duane et al., 1987)) within the Bayesian framework makes it possible to estimate all the unknown parameters (θ j ).The parameters are estimated by updating the prior density distributions P(θ j ) with the observations (D) through the likelihood.This results in a posterior distribution P(θ j |D) for each parameter.The algorithm draws samples from the posterior distribution to provide an average estimation of the parameters with relative uncertainty (Annis et al., 2017;Gelman et al., 1995;McElreath, 2020) Gelman & Rubin, 1992;Vehtari et al., 2021), where 1 ≤ R < 1.05 indicates chain convergence.To ensure efficient algorithm convergence, we chose 1,000 samples, which determine the algorithm's ability to estimate the unknown parameters at the extreme tails of the probability distribution, and four independent Markov chains to sample from.Figure 1 (step 3) shows, as an example, the parameters' prior probabilities, and the posterior probability for one event and four chains, which converge on the same distribution.

Assessment of Model Performance
We here present a general strategy for assessing model performance and evaluating the extracted parameters, which can be applied to any constituent and adapted according to the available data set and analysis goals.To simplify the discussion, we will refer to the parameter vectors introduced in Equations 6-9 as p, h, l, q l .Of these, q l was included in the model because it was required by the mass balance proposed by House and Warwick (1998).In the original work, q l was measured independently and used to find the other parameters (House & Warwick, 1998).However, these measurements were not available in our study: to obtain an unbiased estimate of q l , we therefore let the model define it from a prior probability and available observations.Although the study of q l can provide relevant information on the groundwater status (Duncan, 2019), it is not particularly informative for studying hysteresis behavior and was therefore excluded from further analyses.After modeling, we re-scaled the parameters to the original standard deviation of constituent and discharge to facilitate the model interpretation.
We considered the 2.5th, 25th, 50th, 75th and 97.5th percentile values to study their behavior under uncertainty.Of these, the 50th percentile represents the expectation of the parameters, that is, the region where most of the data is expected to lie, while the other percentiles provide information on the shape and the spread of their probability distribution.Therefore, we considered the 50th percentile values of the parameters as the point estimates (Chakraborti & Li, 2007;McElreath, 2020), for which we evaluated the model performance.

Water Resources Research
10.1029/2023WR035427 We utilized quantitative metrics to evaluate the model's ability to reproduce hysteresis, including its size, rotation and slope.We refer the reader to the Supporting Information S1 for a detailed definition of these metrics.We used the root mean squared error (RMSE) to evaluate the modeling error on individual loops and on the time series.To ensure a scale-invariant comparison among hysteretic loops, we utilized normalized values of RMSE (i.e., NRMSE).Low RMSE and NRMSE values indicate that the model is statistically appropriate to describe constituent concentration during hysteresis.We evaluated the size, rotation and slope of hysteretic loops by comparing well-established hysteresis indices on observed and modeled loops.We used the hysteresis index HI (Lloyd et al., 2016a;Vaughan et al., 2017) to quantify the size and rotation of hysteretic loops and the flushing index FI (Vaughan et al., 2017) to quantify the slope.Both indices are defined in the range [ 1,+1]; a negative flushing index indicates that dilution predominates over enrichment, and a positive HI index indicates clockwise rather than anti-clockwise hysteretic loops.We considered the indices computed on the observed loops as a reference for comparison, discarding from further analysis events where HI could not be calculated because a loop was not delineated.We used Spearman's rank correlation coefficient ρ (i.e., "rho") to measure the similarity of the indices on observed and modeled loops, and RMSE to quantitatively estimate the deviation in loop size and slope.Finally, as an additional metric to assess the model's ability to reproduce the rotation and slope of hysteretic loops, we calculated the percentage of loops with the same sign.
Low RMSE values and high Spearman's rank correlation coefficients indicate that the model accurately reproduces the size, rotation and slope of hysteretic loops.

Statistical and Temporal Analysis of Model Parameters
We performed a statistical and temporal analysis of the model parameters to verify that extracting them from observations allows us to describe the relevant hydrochemical properties established by the mass balance.First, we statistically assessed the ability of the parameter p to describe the size and rotation of hysteresis.We then evaluated the extent to which the model can describe the temporal behavior of diffuse and low-flow inputs.

Size and Rotation of Hysteresis:
We compared p with the hysteresis index HI calculated on observed loops.Note that while the signs of the two indices are directly comparable, their scales are not.The physical limits of p depend on the concentration and discharge values being analyzed and are not necessarily constrained in the range [-1, +1] as for HI.In general, strongly positive or negative values of the parameter p characterize a high sensibility of constituent concentration to water discharge variations (House & Warwick, 1998).We compared the indices with Spearman's rank correlation coefficient, and we calculated the percentage of loops with matching signs.A high percentage and Spearman's rank correlation coefficient indicate a good ability of the parameter p to describe the size and rotation of hysteretic loops.

Diffuse and Low-Flow Contributions:
We compared the parameter h with the flow-weighted concentration (C f , Equation 3), and l with the average constituent load observed in each event to evaluate the contribution of diffuse and low-flow inputs to the corresponding total-flow observations in each event.We quantified these by average percentage contributions, as specified in the Supporting Information S1.It should be noted that the flow percentages associated with diffuse and low-flow inputs are not known in advance from House and Warwick (1998).However, these estimates allow an understanding of how much diffuse and low-flow inputs contribute to the total constituent concentration and load during hysteresis and can be compared with existing knowledge to determine if the model parameters retain their significance.
We evaluated the temporal trends of the extracted parameters to understand the long-term behavior of diffuse and low-flow inputs.We used the non-parametric Hamed and Rao's Mann-Kendall test corrected for autocorrelation (Hamed & Ramachandra Rao, 1998) (significative lag = 1) to statistically test for the presence of a monotonic increasing or decreasing trend.Finally, we correlated the extracted parameters with hydro-climatic variables using principal component analysis (PCA) to verify further that diffuse and low-flow contributions are correctly described and to understand the influence of natural factors on them.

Principal Component Analysis (PCA)
PCA is a multivariate dimensionality reduction technique which allows one to capture the highest variance in the data, extracting the most relevant patterns and correlations between the variables involved in each hysteretic event.We selected relevant hydro-climatic variables, namely water and air temperature (T w(5d) , T air(5d) ), conductivity (at 20°C) (Cond), pH, oxygen and precipitation depth (API7, API14 (Onderka et al., 2012)), together Water Resources Research 10.1029/2023WR035427 with additional features of hysteretic events, such as average values of constituent concentration and load (C (mean) , Load (mean) ) and discharge (Q (mean) ), or the ratio of coefficients of variation (CV ratio (Musolff et al., 2015)).Supporting Information S1 reports details on the considered variables.Only events with p agreeing in sign with HI were considered in the PCA, to evaluate the interaction of hydro-climatic factors on the extracted parameters only for reliable loops.

Chloride as a Case Study: Sampling Site and Data Set
This study uses data stored in the RIWA-Rijn database (Nieuwegein, NL) for the measurement station of Lobith, located on the lower river Rhine, near the Dutch-German border.We tested BaHys on daily time series of water discharge (Q) and chloride concentration (Cl) for a total period of 28 years and 10,227 time points.Table 1 reports detailed information and statistics about the basin and the considered data set.We retrieved the meteorological time series from the Royal Netherlands Meteorological Institute (KNMI) website (KNMI, 2021), for Volkel weather station, located 40 km (air distance) from Lobith.The data employed in this analysis were made available in a public repository (Cairoli & Souza, 2024).We interpolated missing values for chloride using cubic spline for one to three contiguous missing points, and dynamic time warping for more contiguous missing points to retain the time series pattern, according to Caillault et al. (2017).In two temporal windows (26/01/1995-16/03/1995; 28/08/1998-02/10/1998) we could not trust the interpolation due to low data quality: the included hysteretic events were always discarded.

Delineation of Runoff Events
To delineate the runoff events we employed the storm detection algorithm proposed by Dupas et al. (2016), following the implementation of Musolff et al. (2021) (Musolff, 2021b).Before running the algorithm, we smoothed the water discharge time series with Savitzky Golay filter (M.J. E Savitzky, 1964) with window length nine and polynomial order three.We then isolated the events with the following settings: the start of an event was defined in a moving window of 2 days based on a discharge increase of 0.5%.The end of an event was defined when the discharge dropped more than 10% after the discharge peak and decreased less than 0.01% within the moving average window.Only events lasting at least 7 days, and above the threshold of 20% of the minimum total discharge, were included.By choosing these settings, we aimed to isolate all the runoff events, irrespective of their nature, to include all the hydrological variability in the analysis.The minimum length of 7 days was selected to allow the model to fit a sufficient number of points when modeling hysteresis, and rarely do we observe events lasting for a shorter period with data at a daily frequency.It should be noted that the analyst should adjust the choice of these settings depending on the data set being analyzed.

Observed Water Discharge, Chloride Concentration and Load
Figure 2 shows the raw time series of chloride concentration, load (load = Cl * Q) and water discharge.At Lobith, the runoff events that can generate hysteresis are mainly driven by rainfall and snow melts events, which determine the highest peaks on the water discharge time series (Figure 2a).The observed time series show a decrease in chloride concentration and load from the year 1990 (Figure 2b), with a slower decrease in chloride load.We can rely on these observations to evaluate the ability of BaHys to increase our understanding of chloridedischarge interactions that are not directly measured.We first report on the model's performance at reproducing chloride behavior on the time series and during hysteresis, including its ability to reproduce the size, rotation, and slope of the loops; we then describe the performance of each extracted parameter, and report the results of the PCA on parameters and hydro-climatic conditions.

Modeling the Temporal Variation of Chloride During Hysteresis
We delineated 376 events from the time series of water discharge, which were used to fit the Bayesian model.For each event, the Markov Chain algorithm successfully converged to the target distributions for all the model parameters, with an average Gelman-Rubin convergence diagnostic R = 1.00.The hysteresis index HI could be calculated for 358 of the delineated events, which were therefore considered for further analyses.The length of these events ranged from 8 to 62 days, with an average of 22 days.
Figure 3a shows the modeled and observed chloride concentration, and Figure 3b 20 of the modeled loops as an example.The variability in the chloride time series is accurately explained by the Bayesian model (R 2 = 0.921),   1), chloride hysteretic loops did not show high variability in shape, slope and size (Figures 3b and 4a, Figure S1 in Supporting Information S1), and predominantly exhibited clockwise dilution behavior.We evaluate if the model parameter p captures this behavior in the following.

Comparison of Model Parameter p With Hysteresis Index HI
The comparison of the parameter p with the hysteresis index HI is displayed in Figure 4b.The indices agree in the sign for 87.7% of the loops (i.e., 314 over 358 loops), and, despite being defined on different scales, present a high positive correlation (ρ = 0.701), confirming that p compares well with the hysteresis index HI in describing the size and rotation of hysteretic loops.Considering the concordant loops, 86.3% resulted clockwise and 13.7% anticlockwise.Being the flushing index FI negative for the majority of the loops (Figure S1 in Supporting Information S1), it is possible to conclude that chloride concentration at Lobith prevalently shows a clockwise dilution behavior, with the parameter p suggesting little variation in hysteresis events over the considered time series, in agreement with the observations and established hysteresis indices.

Temporal Insights Into Diffuse and Low-Flow Chloride Contributions
Figure 5 shows the temporal behavior of the extracted model parameters h and l with corresponding observations defined for the total discharge.In the long term, the modeled diffuse and low-flow inputs of chloride (Figure 5, Figure S2 in Supporting Information S1) show similar trends and variations compared to their respective totalflow contributions, with the parameter h averaging 62.5% to the total flow-weighted concentration, and the parameter l 51.4% to the total load.Both parameters show a statistically significant decreasing trend over time, with a much stronger decreasing trend for the load at low flow: the trend slope for h equals 0.0565 and the trend slope for l 378.These results suggest that diffuse and low-flow chloride inputs at Lobith contribute similarly to their respective total-flow contributions, with generally higher diffuse inputs during hysteresis, and show a significant reduction of chloride concentration in groundwater and sub-surface runoff over the years, which is also detected, to a lower extent, for diffuse inputs.These relationships are further explored by establishing correlations between the model parameters and antecedent hydro-climatic conditions with the Principal Component Analysis.

Influence of Hydro-Climatic Conditions on Temporal Chloride-Discharge Hysteretic Patterns
The outcome of the PCA analysis is displayed in Figure 6, for 285 of the 314 loops modeled from 1978 to 2003 (precipitation depths were only available from the year 1978).The first three principal components (PCs) capture 57% of the total variability in the system, each of them explaining an independent source of information.The PCA analysis highlights relevant correlations between the model parameters and available observations.The extracted model parameter h is found to follow the same dilution pattern and variation compared to the total flow-weighted chloride concentration (Figures 6a and 6b), statistically confirming what was observed in Figure 5a.A positive For each event, the date of the peak discharge was considered; the moving average computed over five events is shown for each parameter and observation, chosen by visual inspection to visualize temporal trends without excessive variability.The raw time series are shown in Figure S2 of Supporting Information S1.For the sake of completeness, the model estimates of the parameter q l are also reported in Figure S2 of Supporting Information S1.
correlation is also observed for the parameter l and the average total load (Figures 6c and 6d), confirming the similar decreasing trend and variation previously observed (Figure 5b).While both parameters h and l mainly correlate with electrical conductivity, temperature, and precipitation depth, the hysteresis term p does not show relevant correlations with any of these factors.Electrical conductivity is positively correlated with chloride concentration and load due to its ionic nature.Water and air temperature explain most of the variability in the chloride-discharge relationship and determine the dilution behavior of chloride (Figures 6a and 6b).In winter, water discharge is generally higher while chloride concentration is lower, indicating increased dilution in higher discharge volumes.
The effect of precipitation explains 12% of the total variability at Lobith (Figures 6e and 6f).The precipitation indices calculated for the 7 and days preceding the events correlate positively with the chloride diffuse contribution (h), and negatively with chloride load at low flow (l), event duration, CV ratio and water and air temperature.These relationships suggest that the events associated with local precipitation are the shortest of the ) are colored red, variables related to water discharge (q l , Q (max) , Q (mean) ) are colored blue, and variables related to chloride load (l, Load (mean) ) are colored purple.Orange is used for conductivity (Cond), light blue for precipitation depth (API7 and API14), yellow for air and water temperature (T w(5d) , T air(5d) ), and green for pH, O 2 , RL, duration, CV ratio .Note that the signs of scores and loadings of PC1 and PC3 have been reversed for ease of interpretation.
Water Resources Research 10.1029/2023WR035427 CAIROLI ET AL. events modeled on the discharge time series and are most likely to occur at low temperatures.In these events, the diffuse chloride contribution dominates over the low-flow contribution, and chloride presents a chemostatic, rather than chemodynamic, export regime, suggesting little variation in concentration compared to water discharge.

Discussion
4.1.Can BaHys Reproduce C-Q Hysteretic Loops in the Long Term?
We proposed BaHys, a Bayesian modeling framework to study concentration-discharge hysteresis on large temporal scales, and we tested it on chloride over 28 years of monitoring in the river Rhine.The underlying mass balance equation, as presented by House and Warwick (1998), was essential for this.We extended the existing least-squares model to a Bayesian model that allows the simultaneous estimation of all the parameters, with relative uncertainty, solely from observed time series.By evaluating the size, rotation and slope of modeled and observed loops with well-established hysteresis indices, we verified that BaHys is well-formulated to reproduce hysteresis in the long term.Specifically, the comparison of the hysteresis index HI on observed and modeled loops showed a stronger correlation than the flushing index FI.The modeled loops generally presented a smoother profile than the observed loops (Figure 3b, Figure S3 in Supporting Information S1).While this does not considerably affect the estimation of the hysteresis index HI, which averages the differences in constituent concentration on the rising and falling limbs, the FI index is highly sensitive to the concentration at the start and the peak of the event (Liu et al., 2021).Hence, a different modeled concentration at the start (e.g., Figures S3c-S3e in Supporting Information S1) or at the peak of the event (e.g., Figures S3a and S3b in Supporting Information S1) results in large variations in the FI index, but in small variations in the HI index compared to the observed loops.
Among the model parameters, we extracted a hysteresis term (i.e. the parameter p) which was shown to reliably reproduce the size and rotation of hysteresis compared to the well-established hysteresis index HI.Both the p parameter and HI index on modeled loops showed a positive Spearman's rank correlation with the HI index computed on observed loops (HI obs ), with the HI index presenting a stronger correlation than p.Such a stronger correlation is mainly because the two indices are computed differently and defined on different scales.While the HI index is adimensional and defined in the range [ 1, +1] (Lloyd et al., 2016a;Vaughan et al., 2017), the parameter p is a response factor which depends on the analyzed constituent concentration and water discharge (House & Warwick, 1998) and is not necessarily limited within this range.Hence, p parameter values may not always show a strong monotonic relationship with HI obs values, leading to a weaker, but still fairly strong (ρ = 0.701), correlation than HI.The two indices also characterized similar percentages of modeled loops concordant with observed loops, indicating that they agree in describing loops' rotation.Note that in this study we expressed p with its unit measure to show the variation of the hysteresis response to the measured concentration and discharge values.For comparative analyses across different data sets, the scaled value of p as returned by the model (i.e., independent from concentration and discharge scales) could be directly employed, similarly to what proposed by Musolff et al. (2021).

How Well Does BaHys Describe Diffuse and Low-Flow Contributions?
BaHys extracted two parameters, h and l, which describe the diffuse and low-flow constituent contributions in the river.Chloride is a geogenic solute, and its highest concentrations are found in groundwater: during runoff events diffuse sources dominate low-flow contributions due to increased dilution of direct runoff in groundwater (Musolff et al., 2021;Rose et al., 2018).The parameters h and l were shown to correctly describe this mixing behavior for runoff events dominated by local precipitation (Figure 6f).Due to their short duration (Figure 6f), these events most likely originated from isolated precipitation events.At Lobith, however, the water discharge is determined by a much higher hydrological variability and is prevalently dominated by upstream snowmelt events from the Alps (Pinter et al., 2006), which overlap with local precipitation, resulting in events with multiple peaks on the time series.Isolating these overlapping peaks was difficult with relatively low-frequency data such as those available in this study, resulting in events lasting several days or weeks.Therefore, the parameters h and l could reliably describe the mixing behaviors for low-flow and diffuse sources only for well-isolated short events.In these events, chloride showed a prevalent chemostatic behavior, typical of geogenic solutes (Godsey et al., 2019; Water Resources Research 10.1029/2023WR035427 Musolff et al., 2021;Rose et al., 2018), indicating that its concentration is rather stable under widely varying water discharge.
Although mixing behavior is extremely important for analyzing the dynamics of constituents in the river, trend analysis can provide equally important insights into the status of surface and groundwater quality.It was found that h and l show similar decreasing trends and variations compared to the corresponding total-flow contributions .Consistent with the chemostatic nature and dominant dilution behavior of chloride, diffuse inputs were not found to be considerably higher than low flow contributions over time.Considering the uncertainty over the percentiles (Figure 5), these findings suggest that our model correctly characterizes average diffuse and low-flow temporal trends.Comparing the outcome of BaHys with a hydrochemical model and including data at a higher frequency would ultimately ensure that these estimates are physically meaningful, and will be considered for future work.

Model Insights Into Long-Term Chloride-Discharge Dynamics
By considering chloride as a case study, we showed that BaHys allows a reliable interpretation of hysteretic loops by extracting relevant parameters.Analysis of these parameters (i.e., p, h, and l) revealed the dominant mixing mechanisms that underlie chloride-discharge hysteresis.The hysteresis term p showed a clockwise behavior for most of the extracted loops over 28 years of monitoring, with dilution dominating over enrichment during these events.Previous studies on chloride report similar findings: Rose et al. (2018), for example, show a clockwise behavior for chloride in 74% of the runoff events analyzed, and Musolff et al. (2021) register a similar behavior for specific conductance, comprising chloride.Both studies related this behavior to the prevalent geogenic origin of chloride: during runoff events chloride tends to present an initial groundwater concentration source, that progressively mixes with precipitation and soil water sources, resulting in a dilution during hysteresis (Musolff et al., 2021;Rose et al., 2018).The observed mixing behavior was also confirmed by analyzing diffuse and low-flow constituent inputs (i.e.parameters h and l) during short runoff events.
Comparing the extracted parameters with available knowledge and (antecedent) hydro-climatic, biogeochemical and anthropogenic factors can reveal the drivers of constituent export in the long term.In the case of chloride, the long-term decreasing trends of diffuse and low-flow inputs reflect the effect of anthropogenic activities in the Rhine.In the past, chloride concentrations increased considerably due to the intense coal mining activity, which polluted the Rhine water for decades until their cessation in the 80 s (IKSR, 2021).The observed decreasing trends in the extracted diffuse and low-flow inputs since the late 80 s thereby indicate that mitigation measures have effectively limited pollution from mining in both surface and groundwater.These findings suggest that mining activities, whose influence at Lobith is still strong due to the tributaries Lippe and Emscher (IKSR, 2021), should be routinely monitored to keep chloride concentration within the permitted limits.
In our study, none of the included hydro-climatic conditions were strongly correlated with the hysteresis term p from the PCA analysis (Figure 6).This suggests that chloride hysteresis behavior does not vary considerably under different hydro-climatic conditions over 28 years.Comparable behaviors were observed by Musolff et al. (2021) for specific conductance across different watershed and hydro-climatic conditions and hint to similar groundwater sources and mobilization mechanisms over time (Musolff et al., 2021).However, several factors can influence chloride hysteresis patterns at Lobith.For instance, chloride occurs naturally as soluble salt in mineral rocks and is introduced by several anthropogenic sources, such as mining, industrial, and agricultural activities (IKSR, 2021).This variability in sources may influence chloride-discharge hysteresis.A more in-depth analysis of the possible sources and drivers of chloride-discharge hysteresis at Lobith is therefore necessary to draw sensible conclusions and suggest mitigation measures for chloride export based on hysteresis patterns.More comprehensive statistical models, for example, path models (Cairoli et al., 2023;Fernandes et al., 2018), could help relate the spatio-temporal influence of different sources to these patterns.However, this goes beyond the scope of the current study.

Potential and Challenges of Bayesian Modeling of Long-Term C-Q Hysteresis
We have used a simplified mass balance to model and investigate long-term concentration-discharge hysteresis.This simplified description was shown to work well for capturing hysteresis under hydro-climatic variability for chloride.The frequency of our data set mainly allowed the description of average behaviors of diffuse and lowflow inputs, which was sufficient to analyze temporal trends.Even with relatively low-frequency data, the mixing behavior could be reliably reproduced for well-isolated short events.However, longer events were also isolated, for which the initial assumption that diffuse and low-flow sources remain constant may be too approximate, hindering the study of relevant interactions.The inclusion of time-varying rather than event-varying parameters or of higher frequency data should therefore be considered to allow interpretation of diffuse and low-flow mixing behavior along the time series.
The mass balance used by BaHys is generalizable to many constituents and higher temporal resolutions.For example, House and Warwick (1998) applied it to analyze nitrate, ammonium, calcium, phosphorus and silicon data at a two-hourly frequency; Bieroza and Heathwaite (2015) used it to describe hysteresis of phosphorus at a 15-min interval.High-frequency data, hourly or sub-hourly, are becoming increasingly available and provide an opportunity to examine in detail the diffuse and low-flow mixing behavior.Bayesian modeling can extend the mass balance in the long term for higher-frequency data.To demonstrate this, we tested BaHys on 15-min frequency data of spectral absorbance at wavelength of 254 (SAC 254 ), made available and analyzed by Musolff et al. (2021) (Musolff, 2021a).SAC 254 was collected at a much higher frequency and showed enrichment rather than dilution behavior during hysteresis (Figure S4 in Supporting Information S1) (Musolff et al., 2021).Although the model did not achieve R convergence for all the parameter values on the time series, BaHys could correctly reproduce SAC 254 -discharge hysteresis along the time series (Figure S5 in Supporting Information S1), including size, rotation and slope of hysteretic loops (Figures S6 and S7 in Supporting Information S1).As for daily chloride, similar relationships were found between the hysteresis index HI and the parameter p (Figure S6 in Supporting Information S1).In this data, the contribution of the parameter h was dominant over the parameter l (Figure S8 in Supporting Information S1), indicating an increased SAC 254 concentration during runoff events, in line with the findings of Musolff et al. (2021).Hence, BaHys provides a flexible framework that may be generalized to several constituents at varying temporal resolutions.In this framework, the choice of the prior density distribution plays a key role.On the one hand, assuming weakly informative priors adds less specific information about the studied problem, but helps model generalization.On the other hand, assuming more informative priors allows specific knowledge about the parameters to be implemented in the model, facilitating the search for feasible solutions.While weakly informative priors were sufficient to achieve R convergence for daily data, this was complicated by the higher variability of higherfrequency data.In this case, more informative priors should be preferred to constrain the model only to plausible parameter spaces.
The implementation of the mass balance in the Bayesian framework makes the model potentially suitable for predictive purposes, providing decision-makers with prospective indications to mitigate contaminants pollution in the river.Associating these predictions with uncertainty estimates becomes particularly relevant because it enables decision-makers to make meaningful evaluations accounting for the possible risks (Moges et al., 2021).

Conclusions
We presented BaHys, an approach that extends the study of hysteresis in the long term to extract relevant hydrochemical properties of the concentration-discharge relationship.This was possible by combining an existing mass balance with observed constituent and water discharge time series with Bayesian modeling.Bayesian modeling allowed us to accurately model hysteretic loops, extract these properties only from the observations and associate them with uncertainty estimates.By testing our model on daily chloride over 28 years, we could verify that BaHys correctly describes the size and rotation of hysteresis by a hysteresis term that compared well with a well-established hysteresis index.By evaluating the temporal behavior of the diffuse and low-flow inputs to chloride concentration, we could verify that these properties could be correctly extracted according to existing knowledge and observations, although a higher temporal resolution would be beneficial to better capture their mixing behavior in the long term.Comparison of the extracted parameters with existing knowledge and Water Resources Research 10.1029/2023WR035427 observations finally revealed the mixing behavior and the effect of anthropogenic and natural factors on chloride event dynamics and long-term trends.
The employed mass balance and Bayesian model were shown to apply to different constituents and data at varying temporal resolutions.According to the available resolution, however, special care must be taken to correctly apply BaHys and interpret the extracted parameters.On the one hand, when working with low-frequency data, the extracted hydrochemical contributions may not capture relevant interactions for all the events, especially if lasting more than weeks or months.In this case, extending the model to time-varying rather than event-varying parameters could be considered.On the other hand, working with higher-frequency data may require additional measures to deal with the increased variability compared to lower-frequency data.Typically, this need is indicated by model diagnostics (e.g., convergence diagnostic R).Carefully evaluating these diagnostics and adapting the model accordingly, for example, by considering more informative priors, is particularly important to trust the model estimates.To our knowledge, no data set with a higher frequency was available for chloride over such a long period as 28 years.The increasingly available high-frequency water quality data offer an excellent opportunity to test the proposed framework under various data settings.
By utilizing already collected time series, our approach provides a unified framework for long-term analysis of C-Q dynamics with limited data requirements.This study aimed to show the ability of BaHys to model long-term hysteresis dynamics.Further extending the analysis to include several constituents and the influence of relevant drivers would provide additional insights into the main controlling factors of C-Q hysteresis patterns, ultimately proposing measures to mitigate constituent export into the river.With further development, BaHys could be implemented in more complex numerical models to predict river constituent concentrations, explicitly incorporating C-Q hysteresis.As such, our approach has the potential to become a viable monitoring tool to inform decision-makers on future mitigation measures to control and ultimately reduce river water pollution, in line with the key requirements of the EU WFD.

Figure 1 .
Figure 1.Workflow of the approach presented in this study.The equation in step 1 is a simplified form of Equation1, presented in the Methods section, and describes the relationship between constituent concentration and water discharge in a single hysteretic event j, depending on the parameters θ = {p, h, l, q l }.The equation in step 2 is a simplified form of Equation4, which defines the relationship between constituent concentration and water discharge for the entire time series, according to the set of parameters Θ = {θ 1 , …, θ J }, which vary between hysteretic events.A more detailed explanation of parameters and equations is reported in the Methods section of the manuscript.

Figure 2 .
Figure 2. Observed time series of water discharge (a), chloride concentration (b), and chloride load (c).The black line represents a moving average calculated over a sliding window of 90 days.

Figure 3 .
Figure 3.Comparison between modeled and observed chloride time series and C-Q hysteretic loops from 1976 to 2003.(a) Observed (black) and modeled chloride concentration considering 50th percentile (red, RMSE = 16.4 mg/l, NRMSE = 0.283, R 2 = 0.921), 25th and 75th percentile (dark gray shade), 2.5th and 97.5th percentile (light gray shade) values of the model parameters.The 28-year time series of chloride was divided into three subplots to show the temporal behavior over the entire period.The black rectangular symbols indicate the time intervals characterizing the hysteretic loops displayed in panel (b).(b) Observed (black) and modeled chloride-discharge hysteretic loops considering 50th percentile (red), 25th and 75th percentile (dark gray shade), 2.5th and 97.5th percentile (light gray shade) values of the model parameters.The right-pointing triangle represents the loop starting point.Loops' concentration and discharge values were normalized to the range [0,1] to facilitate visual comparison.

Figure 5 .
Figure 5. Modeled diffuse and low-flow chloride contributions.For each contribution, the 50th percentile value (bold line), and the 2.5th, 25th, 75th and 97.5th percentile values computed with Bayesian modeling for each hysteretic event are shown.Panel (a) h (mg/l) contributes to the diffuse chloride inputs and is compared with the flow-weighted chloride concentration (C f , black line).Panel (b) l (kg/s) indicates the load at low flow and is compared with the average observed load in each event (black line).For each event, the date of the peak discharge was considered; the moving average computed over five events is shown for each parameter and observation, chosen by visual inspection to visualize temporal trends without excessive variability.The raw time series are shown in FigureS2of Supporting Information S1.For the sake of completeness, the model estimates of the parameter q l are also reported in FigureS2of Supporting Information S1.

Figure 6 .
Figure 6.Principal component analysis (PCA) of modeled hydrochemical parameters and antecedent hydro-climatic variables for 285 hysteretic events and three principal components.Panels (a), (c), (e) report the PCA score plot for each principal component, Panels (b), (d), (f) the respective loadings.The colored bars in the loadings refer to the most relevant variables in each component: variables related to chloride concentration (h, C f , C (min) , C (mean)) are colored red, variables related to water discharge (q l , Q (max) , Q (mean) ) are colored blue, and variables related to chloride load (l, Load (mean) ) are colored purple.Orange is used for conductivity (Cond), light blue for precipitation depth (API7 and API14), yellow for air and water temperature (T w(5d) , T air(5d) ), and green for pH, O 2 , RL, duration, CV ratio .Note that the signs of scores and loadings of PC1 and PC3 have been reversed for ease of interpretation.

Table 1
Details of Sampling Site and Sampling Data Set Schick et al. (2018)mSchick et al. (2018).CAIROLI ET AL.