Temporal Variability in the Response of a Linear Time-Invariant Catchment System to a Non-Stationary Inflow Concentration Field

Predicting the effects of changes in dissolved input concentration on the variability of discharge concentration at the outlet of the catchment is essential to improve our ability to address the problem of surface water quality. The goal of this study is therefore dedicated to the stochastic quantification of temporal variability of concentration fields in outflow from a catchment system that exhibits linearity and time invariance. A convolution integral is used to determine the output of a linear time-invariant system from knowledge of the input and the transfer function. This work considers that the nonstationary input concentration time series of an inert solute to the catchment system can be characterized completely by the Langevin equation. The closed-form expressions for the variances of inflow and outflow concentrations at the catchment scale are derived using the Fourier–Stieltjes representation approach. The variance is viewed as an index of temporal variability. The closed-form expressions therefore allow to evaluate the impacts of the controlling parameters on the temporal variability of outflow concentration.


Introduction
Agricultural chemicals are frequently mixed into shallow soil layers. Runoff through agricultural land may pick up soil chemicals and introduce them into surface waters. During the periods of heavy rainfall, the soil chemicals in runoff water can pose a potential threat to the quality of the environment and water supplies [1][2][3]. Modeling of the transfer of agricultural chemicals into surface runoff is therefore needed in order to predict the development of the contaminant plume.
The field evidence indicates very irregular patterns of rainfall in time and, in turn, very complex runoff production in a catchment during a hydrologic event [4]. It therefore requires a large quantity of measurements for accurate predictions of the hydrological response of the catchment. This situation is, however, far from that encountered in practice. The amount of information on all coefficients and parameters appearing in the model for accurate predictions is never enough, so one has to accept a level of uncertainty. As such, the runoff and associated solute transport are often formulated as stochastic processes in order to reflect the uncertainty of our knowledge [5][6][7][8][9][10][11]. Since the hydrological variables are regarded as random, predictions are therefore made in terms of probabilities (or moments) rather than in the traditional deterministic framework [12][13][14][15].
Many contamination incidents are geographically isolated and their sources can be considered as point sources of pollution. Although locally important, contamination from these sources may not be Furthermore, this work considers a catchment with a variable input, Ci(t), and an output, C(t), as a linear time-invariant system. If the system represents a convolution, C(t) and Ci(t) may then be connected by the integral form of [26][27][28] C(t) = t 0 g(t − s)C i (s)ds (1) where g(t) is the transfer function (travel time), t-s is interpreted as the time lag between the input and the output time and s is used as a dummy variable of integration, which disappears after the integral is evaluated. The convolution integral operation may be thought of as a filtering operation on the input, where the transfer function is acting as the filter. Equation (1) treats the input flow as uniform over the catchment and ignores the internal spatial variation of catchment flow. The variables Ci(t) and C(t), as considered here, are represented as temporal fluctuations in inflow and outflow concentrations of an inert solute, respectively, and regarded as stochastic (random) processes in time to account for their irregular temporal variation and the uncertainty of their distribution. This study uses the variance as an index of variability of a random variable C(t) in time.
The transfer function g(t) in Equation (1) is a physical measure that integrates flow path heterogeneity and is directly related to internal catchment processes [29]. For simplicity, its temporal distribution is treated as deterministic. That is, it represents an ensemble mean (an average). Following the parameterization proposed by Kirchner et al. [7,8] and Godsey et al. [30], the present study uses the gamma distribution to represent the temporal behavior of g(t).
where α is a shape parameter, τ m is the catchment mean travel time constant, and Γ(α) denotes the incomplete gamma function of argument α. Note that the time required for water to travel through a catchment is a physical measure that integrates flow path heterogeneity. It has become apparent from previous studies [5,31,32] that the convolution integral model based on travel time distributions of solute migration may provide an alternative way of representing the large-scale transport process. As mentioned by Kirchner et al. [7], the gamma travel time distribution with two parameters is both conceptually and mathematically a more suitable representation of the resulting solute mixing processes as its higher degree of freedom offers greater flexibility in accounting for different flow path distributions.
The α in Equation (2) reflects the heterogeneity in the flow path lengths and velocities of catchments [30]. The smaller the value of α, the greater the variability in the catchment flow path length and velocity compared to the mean. For α > 1, the gamma distribution rises to a peak and then falls off, similar to a typical storm hydrograph, which is often used to model rainfall-runoff relationships. For α = 1, the gamma distribution is a close approximation of the exponential distribution, implying that the mixing of a solute plume only occurs at the sampling point; thus, there is no exchange between flow lines. For α < 1, however, the gamma distribution has a completely different shape, having maximum weight at lags near zero and a relatively long tail. These characteristics represent problematic contaminant behavior, with rapid release of an intense contaminant spike followed by persistent lower-level contamination far into the future.
It is furthermore assumed that the temporal distribution of Ci(t) in Equation (1) satisfies the lumped solute transport model [33][34][35][36] as 4 of 13 where τ C represents the solute response time for the inflow solute-transport system and ξ(t) is a zero-mean stationary stochastic forcing (input concentration perturbation) in time. The initial inflow concentration for the lumped solute transport model considered here is a deterministic constant so that the initial condition for the mean-removed solute transport in Equation (3) is As such, the fluctuations in inflow concentration are contributed solely by the random forcing ξ(t). In addition, the temporal variability of the process can be described by an exponential form of covariance function with the associated spectrum S ξξ (ω) defined as [13,33] S ξξ (ω) = Here, ω is the frequency, σ 2 ξ is the variance of the input ξ field and λ is the correlation scale. The correlation scale measures the temporal persistence of ξ field, which is in the units of time. A larger correlation scale indicates longer temporal persistence. The reader may be referred to Mann and Lees [37] and Schulz and Mudelsee [38] for details on identification of the correlation scale of a time series.
Note that the solute transport equation such as Equation (3) is simply a lumped equation of solute mass balance. It describes the time variation of concentration at a catchment scale. The use of a lumped transport model is justified when the whole catchment is treated as a single hydrological unit so that the spatial variability of input, output and parameters are ignored. As such, the solute transport process is treated as a continuous time series rather than a time-spatial process. In fact, in statistical physics Equation (3) is referred to as the Langevin equation [39] which is widely used to describe the time evolution of a particle behavior undergoing a brown motion.
Equation (1) is the starting point for development of the variance of output concentration, along with Equations (2)- (4). Toward this, the stochastic methodology is applied.

Variance of Inflow and Outflow Concentrations
Starting with the use of Fourier-Stieltjes representations [40] for the perturbed fluctuations in Ci(t) and C(t), namely into Equation (1) leads to In Equations (5) and (6)-(8), Λ C i and Λ C are termed as the modulating functions [40] of C i (t) and C(t) processes, respectively, ω denotes the frequency, and Z(ω) represents an orthogonal random process, in the sense that the increments dZ(ω) and dZ(ω ) (at any two distinct points ω and ω ) are uncorrelated random variables, Appl. Sci. 2020, 10, 5356 where * denotes the complex conjugate. By comparison, Equation (8) admits a Λ C − Λ C i relation of the form Furthermore, by virtue of the orthogonality of Z(ω) and the representation theorem for C(t) (namely, Equation (7)), the variance of outflow concentration takes the form where E[ dZ(ω) 2 ] represents the spectral content of fluctuations of the random ξ field in the frequency range (ω,ω + dω). It can be shown that [40] Inferring from the representation theorem for σ 2 ξ , and Equations (11) and (12), it follows that the evolutionary power spectrum for the outflow concentration process S CC (t; ω) can be expressed as That is, upon the solution of Equations (11) and (14) we can evaluate Equation (10) and, therefore, Equation (8).
The solution of Equations (3) and (4) in the Fourier frequency domain can be found by substituting Equation (6) and the Fourier-Stieltjes representation for the stationary ξ process [24,40], into Equations (3) and (4): where ς = 1 + iωτ C and η 1 = exp(−t/τ C ). Within a similar framework for evaluating S CC (t; ω), the evolutionary power spectrum for the inflow concentration process S C i C i (t; ω) and the corresponding variance of inflow concentration σ 2 C i can be computed as follows: Combining Equations (12) and (16), the integrand in Equation (17) can be represented as Appl. Sci. 2020, 10, 5356 6 of 13 The variance of inflow concentration is then found, from substituting Equation (18) into Equation (17) and integrating it over the ω domain, to be where In the limit of t → ∞, η 1 approaches 0 and, therefore, Equations (18) and (19) approach finite values as i.e., the corresponding C i process becomes stationary. It can be inferred from Equations (18) or (19) that the variance of inflow concentration induced by a random ξ persists longer for a larger value of λ or τ C .

Results and Discussion
With Equations (2), (5), (8)-(10), and (16), we are in the position of evaluating the variance of outflow concentration by combining Equations (13) and (14). The evaluation of the integral in Equation (11) cannot be performed analytically for the case of α in Equation (2) being a fractional number. However, to take the advantage of closed-form solutions for Equations (13) and (14), two cases are considered: (1) α = 1 and (2) α = 2. These expressions will be used to investigate the influence of various parameters on outflow concentration variability. It is expected that the concentration variability behavior for the case of α being a fractional number will be qualitatively similar to that of being an integer, although not quantitatively.
For a quantitative evaluation of the analytical developments presented in this section, it is convenient to present the analytical results as the following dimensionless variables: t/λ (dimensionless time), τ m /λ (dimensionless catchment mean travel time), τ C /λ (dimensionless solute response time), τ m /τ C (a ratio of the catchment mean travel time to the solute response time), t/τ C (= (t/λ)/(τ C /λ), a ratio of the dimensionless time to the dimensionless solute response time), and τ C ω (a dimensionless frequency). Note that dimensionless graphs provide much more information than when dimensions are included because it is possible to cover a wider range of the parameters. Figure 1 illustrates the effect of the inflow response time scale τ C on the variance of outflow concentration (Equation (24)) at various time scales. It is seen that the inflow with a larger τ C will lead to higher outflow concentration variability. An increase in τ C enhances the inflow concentration variability and, in turn, the outflow concentration variability. i.e., the temporal variability of outflow concentration persists longer for larger values of τ C . convenient to present the analytical results as the following dimensionless variables: t/λ (dimensionless time), τ m /λ (dimensionless catchment mean travel time), τ C /λ (dimensionless solute response time), τ m /τ C (a ratio of the catchment mean travel time to the solute response time), t/τ C (= (t/λ)/(τ C /λ), a ratio of the dimensionless time to the dimensionless solute response time), and τ C ω (a dimensionless frequency). Note that dimensionless graphs provide much more information than when dimensions are included because it is possible to cover a wider range of the parameters. Figure 1 illustrates the effect of the inflow response time scale τ C on the variance of outflow concentration (Equation (24)) at various time scales. It is seen that the inflow with a larger τ C will lead to higher outflow concentration variability. An increase in τ C enhances the inflow concentration variability and, in turn, the outflow concentration variability. i.e., the temporal variability of outflow concentration persists longer for larger values of τ C .  The variance of outflow concentration represented as function of τ m is depicted in Figure 2 using Equation (24) for various values of λ/τ C . The parameter τ m apparently plays the role in smoothing out the concentration variability. This can be attributed to the reduction in the persistence of outflow concentration perturbation with τ m which reduces the concentration variability. Here, the persistence is referred to the tendency for low values to be followed by low values and high values by high values. Additionally, from Figure 2, we see that for a fixed τ m , an increase in λ results in higher variability of outflow concentration. A larger value of λ increases the persistence of inflow concentration and, in turn, the outflow concentration, leading to higher concentration variability.   For α = 2, the transfer function g(t) in Equation (2) has the form which leads Equation (14) eventually to where µ 4 = 4 + τ 2 m ω 2 and expressions for Γ 1 − Γ 6 are given in the Appendix B. It follows from Equation (25) that the associated variance of outflow concentration is obtained as where expressions for Φ 1 − Φ 8 are given in the Appendix C. Figure 3 shows a comparison between the outflow concentration evolutionary power spectra given by Equations (23) and (26) at a given time. Evolutionary power spectra have the interpretation of decomposition of total energy (or the strength of the variations) over frequency at time t [40]. It is apparent from Figure 3 that high-frequency variations in the temporal input series are attenuated. A similar behavior is shown in Kirchner et al. [7,8] for the case of stationary outflow concentration time series. The figure also implies that the persistence of the concentration field is enhanced by a smaller value ofα. A smaller α will introduce greater variability in outflow concentration.  The parameter α in Equation (2) could be a fractional number. In this case, the integral in Equations (11) and (13) may not be evaluated analytically. Since the transfer function is acting as a filter, it is expected that the outflow concentration variability behavior for the case of α being an integer will be qualitatively similar to that for the case of α being a fractional number, but not quantitatively.
The determination of the transfer function in Equation (1) is of great importance in an application of the convolution integral model for modelling of hydrological processes. A brief introduction of using the cross-correlation technique in determining the transfer function is given as follows. Using Equations (6) and (7), it follows that The parameter α in Equation (2) could be a fractional number. In this case, the integral in Equations (11) and (13) may not be evaluated analytically. Since the transfer function is acting as a filter, it is expected that the outflow concentration variability behavior for the case of α being an integer will be qualitatively similar to that for the case of α being a fractional number, but not quantitatively.
The determination of the transfer function in Equation (1) is of great importance in an application of the convolution integral model for modelling of hydrological processes. A brief introduction of using the cross-correlation technique in determining the transfer function is given as follows. Using Equations (6) and (7), it follows that (28) where R CC i the cross-correlation between C i (t) and C(t) processes. Following Priestley [40], R CC i has a Fourier-Stieltjes representation of the form where S CC i is the cross-spectrum between C i (t) and C(t) processes. A comparison of Equation (28) with Equation (29) yields or At large times, Equation (31) reduces to where Given finite observations on the input (inflow concentration) and output (outflow concentration) processes, Equation (32) enables us to estimate the form of the transfer function by replacing the terms in the numerator and denominator by their respective sample estimates.
Most of the stochastic literature on the analysis of output time series from a linear time-invariant catchment system treats the input to the system as a stationary process, while this work is concerned with a temporal nonstationary input process, which leads to a nonstationary output process from the catchment system. In addition, the analyses presented in the above-mentioned references are carried out using the field observations (or experimental data) as an input time series to the system. In this work, the input nonstationary process is generated from solving a lumped solute transport model (the Langevin equation). That is, the present analysis is carried out using a combination of the convolution integral and Langevin models. To our knowledge, the closed-form expression for the variance of outflow concentration from the catchment system, such as Equation (24) or Equation (27), has not been presented in the literature so far using a nonstationary process as an input to Equation (1). The expression provides a useful way to evaluate the influence of parameters appearing in the models on the variability in outflow concentration. It is also interesting to note that the closed-form expression for the variance of nonstationary process with respect to a lumped solute transport model given by Equation (19) has never been presented in the literature before.
The proposed approach has the advantage of operating at the catchment scale and therefore, if an appropriate input-output transfer function can be specified at that scale, predicting the response of the catchment to any inflow concentration does not require specification of the small-scale spatial heterogeneity of precipitation and soil properties.

Concluding Remarks
The stochastic modeling approach recognizes that hydrological variables are affected by heterogeneity and regards them as random. This randomness leads to defining models of flow and transport in a stochastic context. An important feature of the stochastic approach presented here is its providing a rational framework for quantifying the variability in outflow concentration from a catchment.
The existing stochastic models for the analysis of solute transport in catchments generally assume stationarity in the concentration fields of inflow and outflow. The stochastic nature of non-stationary rainfall-runoff processes, and thus the solute transport processes through the runoff are ignored. The results of the present research, which deals with the analysis of the non-stationary catchment solute transport process, should be valuable for a better understanding of the behavior of the transport of pollutants in surface runoff at the outlet of natural catchments. Import assumptions of this study are distributed contaminant sources and a linear rainfall-runoff process. These assumptions allow to neglect the spatial variability of the hydrological properties of the catchment and to express the output concentration of the catchment as a convolution of the input concentration and the transfer function.
In this study, the closed-form expressions for the variances of inflow and outflow concentrations of inert solutes at catchment scale are developed using the Fourier-Stieltjes representation approach. Based on the closed-form solutions, the impacts of the controlling parameters on the temporal variability of concentration are investigated. The results of the analysis indicate that the solute response time for the inflow solute-transport system and correlation scale of input concentration field play the role in increasing the persistence of outflow concentration perturbation, while the catchment mean travel time scale reduces it. A higher persistence means a less fluctuation in the outflow concentration around the mean, which leads to a larger concentration variability. In addition, the shape parameter reduces the persistence of outflow concentration perturbation and, in turn, the outflow concentration variability.