Using the backward probability method in contaminant source identification with a finite-duration source loading in a river

Violation of industries in discharging their effluents into rivers leads to river pollution, which endangers the environment and human health. Appropriate tools are needed to deal with violations and protect rivers. The backward probability method (BPM) is one of the most recommended tools identifying the release time and location of the pollutant source. However, the BPM generally was developed for groundwater and spill injection. Since most industries inject their effluents with a constant rate for a finite duration, the use of prevailing models will have some errors. In this study, a numerical model was developed that could simulate a source with either a finite-duration or spill injection. This model is verified for two hypothetical cases and one real case. The results show that the model can accurately identify the release time and location of the pollutant source.


Introduction
Water is one of the most important and fundamental factors for living organisms. Surface water resources, because of their availability and easy extraction, are widely used in industry, agriculture, and municipal services. Surface water plays a prominent role in human societies (Cheng and Jia 2010;Jing et al. 2020;Lu et al. 2015;Mazaheri et al. 2015;Rasheed et al. 2019;Saha et al. 2017). Increasing the rate of using surface water resources led to increasing pollution of these resources and scarcity of them (Amiri et al. 2019). Scarcity in water resources affects human lives and threatens economic and social development (Duan et al. 2016;Jing et al. 2020). Therefore, protecting surface water resources has become a global and security issue (Mekonnen and Hoekstra 2018).
Environmental contaminant transport problems are broadly categorized into two subgroups, source-based, and receptor-based problems. In source-based problems, information about the source is available, and the goal is to determine the concentration downgradient from the source (Neupauer and Wilson 1999). In order to respond to these problems, researchers use the solution of the transport model (forward model) in which the advectiondispersion equation (ADE) is the base for analysis (FISCHER et al. 1979;Runkel 1996;Van Genuchten 1981). In receptor-based problems, contaminants are observed at monitoring points. In these cases, the main point is identifying source location, release time, and injected mass. Researchers (Boano et al. 2007;Ghane et al. 2016;Neupauer and Wilson 1999;Tong and Deng 2015) have used the inverse methods in order to respond to these problems. The approaches that are used in the inverse methods can be classified into three categories: simulationoptimization approach, mathematical approach, and probabilistic approach (Mazaheri et al. 2015).
The simulation-optimization approaches use a combination of optimization algorithms and the forward model to solve the inverse problem. This method, despite its simple formulation, requires high computational power. Alapati and Kabala (2000) used the least-squares method without any regularization methods in order to reconstruct the released mass in groundwater aquifers. Mahinthakumar and Sayeed (2005) claimed that combination methods have a better performance than basic optimization, in which the goal is finding the relative minimum. They combined the genetic algorithm with the basic optimizations in order to identify the source location. Xing et al. (2019) proposed an ensemble surrogate model to identifying groundwater contaminant sources. This model was following three individual surrogate models: kriging, radial basis function, and Monte Carlo method. They could effectively identify not only conservative contaminant source but also contaminants containing chemical reaction.
Mathematical approaches solve inverse problems in a direct manner and reduce computational costs. However, complexity is considered a disadvantage of these approaches. Skaggs and Kabala (1994) applied a mathematical method to calculate the release rate of pollutant sources; actually, they solved the Fredholm integral equation. Thus, they used the Tikhonov regularization technique in a one-dimensional domain to solve the inverse form of the ADE for groundwater. Milnes and Perrochet (2007) used the concepts of transfer functions theory to find the location of pollutant sources in groundwater aquifers. They could simulate this phenomenon as an inverse problem by using concentration at the detection point. They showed that in the inverse solution, concentration contours start to gather until these contours become a point that represents source location. Mazaheri et al. (2015) tried to determine the release rate of pollutant sources in the river with non-steady non-uniform flow conditions. They tried to resolve the problem of ill-posed equations by using the least square, the Green's function, and the Tikhonov regularization method.
Probabilistic approaches use probabilistic distributions, and one of the important advantages of this approach is small simulations and computations (Amiri et al. 2019). Liu (1995) established a continuous inverse model that could identify the location of spill-releasing pollutants into groundwater aquifers. In order to produce a probability density function (PDF), he used Laplace transformations and inversed the flow velocity. Neupauer and Wilson (1999) tried to use the adjoint method to solve the inverse ADE and produced a PDF. They verified their model using the model that Liu (1995) presented. Neupauer and Wilson (2001) enhanced their method and showed that the adjoint method could be used for multidimensional problems. They showed how the numerical methods could be combined with the adjoint method in order to identify the release time and source (Neupauer and Wilson 2004). Cheng and Jia (2010) were one of the first to try to identify the release time and the source location in a river. They used the adjoint method under steady non-uniform flow conditions; also, they assumed that the pollutant source is an instantaneous point source of a contaminant. Ghane et al. (2016) realized the release time and the source location in a river using the adjoint method by assuming one dimensional, steady, nonuniform flow conditions, and an instantaneous point of source contaminant. They used a hypothetical example and real river conditions in order to verify their backward model. Jing et al. (2020) tried to obtain the pollutant release mass in a river by using the adjoint method. They assumed that the pollutant source is a finite-duration source, not instantaneous; their model was established for one-dimensional, steady non-uniform flow conditions. They verified their model using two hypothetical examples and a real case study.
Rivers are known as the main source of water for industry, irrigation, and domestic use. Therefore, they must be constantly monitored for quality; in this case, identifying the source of the pollutant in the rivers is important. Traditional pollutant source identification methods require a lot of data monitoring, including water quality, topographic, and hydrological data, while numerical models of water quality modeling can be more efficient. The backward probability method as a numerical model can identify contaminant sources with limited data. In previous researches, the backward probability method is used for source identification in groundwater with an instantaneous point source. It is clear that many industries discharge their pollutants into the river in a step-loading manner. Therefore, to identify the industries that violate, it is necessary to develop a numerical model appropriate to this type of injection. Furthermore, the instantaneous and finite-duration injection should be discriminated. In this study, the adjoint method is applied for a source with a finite-duration injection, and the adjoint equations for steady uniform flow conditions and conservative pollutants are derived. The proposed model is validated by hypothetical and experimental examples.

Study area
In order to verify the proposed model, a hypothetical problem and a real case are considered.

Hypothetical area
In the hypothetical case, a hydrodynamically acceptable example should be defined. In this case, river width, flow velocity, roughness coefficient, and river slope were assumed; with these assumptions, the water depth and dispersion coefficient were calculated using Manning's equation and Fisher's equation (Hessel et al. 2003;Fisher, 1975), respectively. Suppose that there is a factory upstream of a straight river. This factory, which is located 21500 m upstream of the monitoring point, discharged its sewage continuously in the river for a finite duration. Fig.  1 shows the factory and the monitoring point location in the river. Table 1 lists the flow and transport properties of this example.

Chicago waterway area
Jackson and Lageman (2013) did a field trace experiment in which Rhodamine WT (RWT) dye was injected for a finite-duration source estimation to the Chicago Sanitary and Ship Canal (Fig. 2). In this canal, six stations measured the dye concentration, two of which, fluorometer FL291 located 9930 m downstream of the source and FL293 located 7033 m downstream of the source, are used in this study. The US Geological Survey collected the data. Hydrodynamic parameters used in the proposed model are the calibration results of Zhu et al. (2017). More information can be found in Jackson and Lageman (2013) and Zhu et al. (2017). Table 2 summarizes the field tracer experiment information.

Methodology and methods
In this section, the forward model and the backward model are described. The prerequisite of solving the backward model is the forward model results. It was stated that the ADE was the base equation in the forward model. In order to solve the ADE, hydrodynamic conditions must be known; thus, a hydrodynamic model is used. Series of inputs and outputs derived from the hydrodynamic model are used in the ADE. The forward model results are used as inputs to the backward model to find the release time and source location. Furthermore, the identification process, the numerical model, and the way that the proposed model was verified are described.

Forward model
The basic governing equation of contaminant transport is the ADE that is shown as Eq. (1). This equation is for one dimensional, steady, and uniform flow conditions (Chapra 1997): where C is the concentration, D is the longitudinal dispersion coefficient, U is the flow velocity, S is the source term, K is the first-order decay rate coefficient, t is the forward time, and x is the longitudinal coordinate. The source term defines the characteristics of the source and shows how pollutants discharge into the river. In this study, identifying the source of the pollutants in the river with a finite duration of injection (step loading, Fig. 3) is the main idea. Therefore, the source term can be defined, as shown in Eq. (2): where C 0 is the initial concentration at the source location, δ is the delta Dirac function, H is the Heaviside function, t 1 is the time that pollutant injection starts, and t 2 is the time that pollutant injection stops. Considering a conservative solute (K=0), the ADE and its initial and boundary conditions became as shown in Eq.
The result of solving equation Eq. (3) is the concentration at time t and location x. Normalizing concentration by the total concentration entering the system defines the forward location probability (FLP) from the source to the detection point (Neupauer and Wilson 1999). Equation (5) shows this expression: Fig. 2 Map of Chicago waterway area Fig. 3 Step loading with a finite duration In this study, it is assumed to have a steady and uniform flow condition. Manning's equation and Fisher's equation are used in order to determine the hydrodynamic parameters (A, U, and D). Equation (6) represents Manning's equation, and Eq. (7) represents Fisher's equation.
where n is a Manning's roughness coefficient, R is the channel hydraulic radius, S is the channel slope, B is the average channel width, H is the average depth, and U * is the shear velocity.

Backward model
The adjoint method is used in order to form the backward model. The adjoint of the ADE for the finite-duration of injection was developed by pursuing the steps that (Neupauer and Wilson 1999) established. The performance measure, P, expresses the state of the system and is defined as Eq. (9): where h(α, C) is a function of the state of the system, and α is a vector of system parameters (U, D, C 0 ). The marginal sensitivity is obtained by differentiating the performance measure with respect to one of the system parameters. Equation (10) shows marginal sensitivity in which differentiation is with respect to C 0 .
where dP dC 0 is defined as a marginal sensitivity, ψ ¼ ∂C is a state sensitivity that indicates the measure of the change in the state of the system. The state sensitivity in Eq. (10) is unknown; thus, in order to calculate the marginal sensitivity, ψ from Eq. (10) must be eliminated by using the adjoint method. The first step in utilizing the adjoint method is differentiating the ADE, and its initial and boundary conditions with respect to the C 0 . Differentiating Eq. (3) with respect to the C 0 gives The next step is obtaining to an equation similar to the ADE in which ψ * , an arbitrary function, is the main variable. It will be shown that ψ * is the adjoint state. Therefore, the inner product of the adjoint state and each term in Eq. (11) is taken after defining the inner product of two functions.
where T is the final time (detection time), x 1 and x 2 are the downstream and the upstream boundaries, respectively. A similar form of the ADE for ψ * is obtained by manipulating each term in Eq. (12).
Substituting Eq. (13), Eq. (14), Eq. (15), and Eq. (16) into Eq. (12) gives The left side of Eq. (17) is equal to zero; thus, it can be added to the right side of Eq. (10) without a change in results. Equation (18) shows the new expression of the marginal sensitivity. It is obtained by rearranging terms after adding the left side of Eq. (17) to the right side of Eq. (10).
The current goal is to eliminate ψ from Eq. (18).
Therefore, it can be assumed that ∂h α;C ð Þ = ∂C 0 ¼ 0, because the system parameters (e.g., U and D) are independent of C 0 . The coefficient of ψ in the second term of integration in Eq. (18) is equal to zero (the adjoint equation). The third and fourth terms are divergence terms and represent initial and boundary conditions for the adjoint equation after integration. These assumptions are shown in Eq. (19), Eq. (20), and Eq. (21).
Substituting conditions that represented in Eq. (9) into Eq. (17) gives As previously indicated, the adjoint state, ψ * , is an arbitrary function and can be defined to eliminate ψ. Therefore, the initial and boundary conditions for the adjoint equation are shown in Eq. (22) and Eq. (23): The final time in the forward model, T, is determined arbitrarily, and t ≤ t d is essential; thus, the detection time can be determined as the final time. In the backward model, backward time is used, τ, which is defined as τ = T − t = t d − t. The definition of the function of the state, h(α, C), is related to the performance measure. When the performance measure was considered as a concentration at a known time and location, i.e., the detection point, the function of the state should be defined in a way that the integration over the time and space domain results in the measured concentration at the detection point. Therefore, the function of the state can be defined as Eq. (24). h α;C ð Þ = ∂C is a Frechet derivative (Saaty, 1981) of h(α, C) with respect to the C. Take the derivative of both sides of Eq. (24) gives Eq. (25). Equation (26) represents the marginal sensitivity. Equation (27) shows the adjoint equation and its initial and boundary conditions.
Results of Eq. (27) are the adjoint state at every point in the time and space domain. Using these results in Eq. (26) gives the marginal sensitivity. Eq. (28) shows backward location probability (BLP).
Identification process The FLP demonstrates the probability of the existence of a contaminant that has been released from the source at a particular location (detection point) after passing the specific time. However, the BLP determines the prior location of a contaminant that has been observed at the detection point (Neupauer and Wilson 2001). Based on the backward probability method, BLP and FLP density functions are equal only at the source location (Wagner et al. 2015). The root mean square error (RMSE) have been used to compare the BLP and FLP density functions (Cheng and Jia 2010; Ghane et al. 2016), which Eq. (29) shows where n is the number of data entry points to be compared, the source location will be identified when the RMSE is minimized. Fig. 4 shows the overall identification process.

Numerical modeling
A numerical model has been developed to compute the forward and backward probabilities. This numerical model is based on the finite volume method with the upwind and fully implicit scheme in space and time. The applied numerical model can solve Eq.

Results and discussion
The proposed model is verified for a hypothetical case and a real case. Also, it is expected that this model could simulate an instantaneous source since, in the step loading source, when the difference between the time that source starts and stops to inject contaminant becomes zero, t 1t 2 =0, the source becomes an instantaneous source. Therefore, this ability was verified by using a hypothetical case.

Finite-duration source (case 1)
After injection, the monitoring point records concentration over time (Fig. 5a), and the numerical model uses Eq. (5) to produce the FLP. The proposed model calculates the BLP for every possible source location and compares it with FLP (Fig. 4b). However, the minimum RMSE between BLP and FLP shows source location. In this case, the model determined x = 30000 m as a contaminant source location since it has the minimum RMSE = 1.89 * 10 −5 . Fig. 5a shows that 6.475 h after injection, peak concentration reaches to the monitoring point. The peak concentration was considered as the presence of a contaminant at a specific location (Neupauer and Wilson 1999). Fig. 5b shows the best match between BLP and FLP, and the highest probability for the release time from the detected source location is 6.475 h before observing contaminants in the monitoring point. As it is clear, the model could identify the release time and location with the least error.
Field tracer experiments (case 2) Fig. 6a shows a comparison between measured concentrations at two stations and the numerical results. As shown, the peak arrival time in both numerical and measured concentrations for the FL291 station is 24.3 hours after injection.

Instantaneous source (case 3)
Case 3 is defined to show that the proposed model can identify the release time and location of the contaminant source when the injection happens as an instantaneous source. In case 3, the hypothetical example that is represented in case 1 is used. However, if the model is used for simulating an instantaneous source, the difference between the time that the source starts and stops injection must be zero, t 2 − t 1 = 0. Fig. 7a shows the concentration time series at the monitoring point; based on the results in Fig. 7a, the FLP became accessible. The peak concentration reaches the monitoring point 5.96 h after injection. Fig. 7b shows the comparison between the FLP and BLP at the source location, which has the minimum RMSE (1.22 × 10 −5 ). Based on Fig. 7b, the highest probability for release time is 5.96 h before the injection. The proposed model has been validated by the three examples. The results show that this model can properly identify source location and release time. Most of the previous studies simulate instantaneous injection. With this in mind, in the third example, the proposed model could simulate instantaneous injection and determine the release time and source location. Furthermore, a sensitivity analysis for hydrodynamic parameters was done, and the results indicated that the model could identify the release time and source location correctly, and changing hydrodynamic parameters does not significantly affect the results.

Conclusions
In this study, an adjoint method was developed for a finite-duration injection to a river in order to expand previous studies. The proposed model successfully simulated steady uniform flow conditions and conservative pollutants. The results showed that when conservative pollutants exist, a monitoring point is sufficient for determining the release time and source location. Furthermore, the model performance in a real case was examined, and it was concluded that the model works accurately. The maximum error for release time and source location, in some cases, was 6%. Also, this model can be used for largescale rivers and instantaneous sources. The adjoint method identifies the source location and release time by simulating the backward model once which is a beneficial ability compared to other methods in the literature. Availability of data and materials There is no availability for data and materials.
Author contribution Hossein Khoshgou is an M.Sc. student, and this paper originated from his thesis, and Seyed Ali Akbar Salehi Neyshabouri is his supervisor and corresponding author.

Declarations
Ethics approval This manuscript is our original, unpublished paper, which we did not submit to another journal before.
Consent to participate and publish It is our great pleasure to submit to you my manuscript entitled "Using the Backward Probability Method in contaminant source identification with a finite-duration source loading in the river" to be published in your admired journal.