ISSN 2306–6172

Volume 5, Issue 1 (2017 ) 53 – 73

NUMERICAL SOLUTION OF AN INVERSE PROBLEM OF DETERMINING THE PARAMETERS OF A SOURCE

OF GROUNDWATER POLLUTION

Karchevsky A.L., Turganbayev Y.M., Rakhmetullina S.J., Beldeubayeva Zh.T.

Abstract The article deals with an inverse problem of determining parameters of ground- water pollution sources. We test three ways to solve the problem on the simulated data for a simple case of contamination. We discover that, in the presence of noise in the data of the inverse problem, the first method does not produce satisfactory recovery results, while the second and third ones are comparable in accuracy of the recovery of required parameters.

Taking into account the ease of implementation, the speed of computing and parallelization feasibility the second method of solving the inverse problem is found to be most preferable.

We also propose a method of finding pollution parameters in general case.

Key words: inverse problem, groundwater pollution, Volterra integral equation of the 1-st kind

AMS Mathematics Subject Classification: 35K20, 35R30

### 1 Introduction

Groundwater plays a highly important role in the water cycle in nature and maintaining the balance of ecosystems. It is a source of clean drinking water for many communities [12]. 75% of the Earth’s surface is covered by water, of which 97% of the world’s water is salty. Of the remaining three percent, two are presented as glaciers, so for the human needs only 1% of world reserves of water is suitable. This amount of water is in the form of surface water and groundwater, which accounts for 96% of the reserves.

General problem of environment pollution is closely related to the problem of groundwater contamination because ground water interacts closely to atmosphere, ground surface and surface water sources [13]. Purification of groundwater pollution is very expensive, so the main method of protection is prevention, that is, identifica- tion of pollution sources and their power, as well as the development of measures to prevent the contaminants getting into groundwater. This applies in particular to the pollution of groundwater by industry and urban settlements that have more or less local character [13, 14]. Contamination of groundwater can occur both from point and distributed sources. In order to determine sources, forms and dynamics of groundwater contamination it requires using of mathematical modeling. There are two approaches to the mathematical modeling of groundwater dynamics: direct and inverse. The direct approach assumes that by solving master equations we calculate unknown parameter such as speed, pressure, impurity concentration, etc. after the initial moment of time at

parabolic. In this regard, it is necessary to highlight the work [31], whose authors have made an outstanding contribution to the theory and numerical methods for the solution of parabolic inverse problems.

There is large quantity of data on the flow, structure and pollution of groundwater by various industrial wastes in some areas of the Republic of Kazakhstan. Analysis of the data and forecast on its basis is usually conducted by methods of information technology (see., for example, [1]-[11]). In this paper, we provide a simple method of data processing, based on the theory and numerical solution of inverse problems. The paper proposes three methods for the solution of the inverse problem under study. We test the methods using simulated data and chose one of them, based on the criteria of the ease of implementation, computing speed and feasibility of parallelization.

### 2 Problem statement

The propagation of groundwater contamination is described by the following equation:

RCt=DCxx−νCx, (1)

whereC(x, t) = ˜C(x, t)−C0,C˜ is concentration of a substance entered into the ground-
water, C_{0} is an average background value of the concentration of the substance, D is
diffusion coefficient, R is retardation factor (R≥1),ν is the velocity of groundwater.

We assume that at the initial time, the concentration of contaminant in groundwater
is equal to C_{0}, i.e.

C(x,0) = 0. (2)

We assume that atx= 0 the impact on groundwater is described by the following boundary condition

−(DC_{x}−νC)|_{x=0} =νf(t). (3)
We assume that at a large enough distance from the point of observation of the
concentration of the substance does not change

C_{x} → 0 (x→ ∞). (4)

For a known functionf(t)and the values D,R and ν it is not difficult to calculate a function C(x, t)that is a solution of direct problem (1)-(4).

In case the function f(t) is unknown we will assume that we have accessible mea-
surements of the concentration of pollution in some remote point x=x_{0}:

C(x_{0}, t) =g(t). (5)

Thus, we can formulate an inverse problem: according to measurements (5), deter- mine the unknown function f(t) if it is known that the reduced concentration C(x, t) satisfies the conditions of direct problem (1)-(4) and the parameters D, R and ν are known.

### 3 Equivalent statements, existence and uniqueness of the solution of the inverse problem

We apply the Laplace transformation to the variable t. Then equation (1) transforms into:

Cˆ_{xx} − ν

DCˆ_{x}− R

DpCˆ = 0, (6)

where C(p) =ˆ L[C(t)] and p is a parameter of the Laplace transformation. Boundary conditions take the form

−(DCˆ_{x}−νC)|ˆ _{x=0} =νf(p),ˆ Cˆ_{x} →0 (x→ ∞). (7)
Additional information takes the following form:

C(xˆ _{0}, p) = ˆg(p). (8)

We find the solution of equation (6) in general form
C(x, p) =ˆ A_{1}e^{λ}^{1}^{x}+A_{2}e^{λ}^{2}^{x}, where λ_{1,2} = 1

2

"

ν D ±

rν^{2}

D^{2} + 4R
D

# .

According to the second boundary condition (7) A_{1} = 0. (8) leads to
C(x, p) = ˆˆ g(p)e^{λ}^{2}^{(x−x}^{0}^{)}.

Substituting this expression into the first boundary condition (7) we obtain the rela- tionship between the Laplace images of functions f(t)and g(t):

ˆ

g(p) = 2ae^{ab} e^{−b}

√

a^{2}+p

a+p
a^{2}+p

f(p),ˆ where a= ν

√RD, b=x_{0}
rR

D. (9)

The inverse Laplace transformation [32, p. 221] for equation (9) leads to the rela- tionship

g(t) =

t

Z

0

K(t−s)f(s)ds, (10)

where

K(t) = 2ae^{ab}
1

√πte^{−}^{b}

2

4t −ae^{ab+a}^{2}^{t}erfc
b

2√

t +a√ t

e^{−a}^{2}^{t}.

We find that the solution of inverse problem (1)-(5) is equivalent to the solution of the Volterra integral equation of the 1st kind (10).

solution of inverse problem (1)-(5) leads to its uniqueness. To the best of our knowledge, the theorem of the existence of solutions of integral equation (10) in general form is still not proven.

However, for practical applications, the contamination could well be approximated by a function [14]

f(t) =βe^{−αt}, (11)

where the parameters α (1/d) andβ (m^{3}/d) characterize the intensity of groundwater
contamination at x= 0.

Suppose the functiong(t)in the interval [0, T]satisfies conditions:

g^{0}(t) +αg(t) = βK(t), g(0) = 0, (12)
where α and β are certain constants. Then the solution of inverse problem (1)-(5) in
the interval [0, T] exists and is given by (11).

Indeed, the solution of Cauchy problem (12) can be written in the form

g(t) =

t

Z

0

K(s)βe^{−α(t−s)}ds,

that coincides with (10) after the change of variables τ =t−s.

It is obvious that the problem of determining of the constants α and β can be formulated as follows: Find the constants α and β, if the function g(t) is known, and we know that it satisfies to the conditions of problem (12).

### 4 Numerical solution of the inverse problem

Showed below variants of the numerical solution of the inverse problem have been tried on simulated data, i.e. some values of the constants α and β that are possible in practice were chosen. Then using function f(t)of form (11) the functiong(t)had been computed, and then a multiplicative random error term was added to it.

g_{δ}(t_{k}) =g(t_{k})

1 +ξ_{k} P
100

,

where ξk is a random variable normally distributed on the interval [−1,1], and P is percentage error value.

Since real data has an error the data smoothing procedure is needed. Based on the properties of the smoothing procedure (good removal of a high-frequency components [37]), we chosed the following two methods for the numerical experiments:

1) moving average

˜

g(t_{i}) = 1
2N+ 1

N

X

l=−N

g_{δ}(t_{i+l})
2) smoothing with a Gaussian kernel:

˜

g(t_{i}) = 1
2N+ 1

N

X

l=−N

G

g_{δ}(t_{i})−g_{δ}(t_{i+l})
b

g_{δ}(t_{i+l})

^{N}
X

l=−N

G

g_{δ}(t_{i})−g_{δ}(t_{i+l})
b

, where

G(t) = 1 0.37√

2π e^{−}

t2 2(0.37)2

and b is the parameter determined be the width of the smoothing window. In this paper, we chose N = 5, b = (2N + 1)τ, where τ = 1 day is the time step, with which we know the function g(t).

We ran a large number of numerical experiments. The results and conclusions were made based on the following pair of test exact solution of the inverse problem:

α= 3.0·10^{−2}, β = 5.0.

Data error of the inverse problem was P = 3%. D = 10(m^{2}/d),R = 1, ν = 1 (m/d),
x0 = 100(m).

In order to calculate the function K(t) it is necessary to calculate the function erfc(t). To do this, we used the expression [38]:

erfc(t) = e^{−t}^{2}(t), (13)

where

(t) = τe^{η(τ)}, τ = 1/(1 + 0.5t),

η(τ) = −1.26551223 + 1.00002368τ+ 0.37409196τ^{2}+ 0.09678418τ^{3}

−0.18628806τ^{4}+ 0.27886807τ^{5}−1.13520398τ^{6} + 1.48851587τ^{7}

−0.82215223τ^{8}+ 0.17087277τ^{9}.

In view of (13) the expression for K(t)can be rewritten:

K(t) = 2a·e^{−}

b

2√ t−a√

t2

√1

πt−a b

2√

t +a√ t

.

### 4.1 The first method of numerical solving

Since the data g(t)must satisfy (12), we offer the following algorithm:

• Step 1: Select the interval[τ1, τ2](see. Fig. 1) from which we take data to recover the parameters α and β;

Figure 1: Function g(t).

• Step 2: Select the time moments t_{1,i} < t_{2,i} (i= 0, N) from the interval [τ_{1}, τ_{2}];

• Step 3: Solve the system of equations

˜

g_{δ}(t_{1,i})α−K(t_{1,i})β = − 1
2τ

˜g_{δ}(t_{1,i}+τ)−g˜_{δ}(t_{1,i}−τ)
,

(14)

˜

g_{δ}(t_{2,i})α−K(t_{2,i})β = − 1
2τ

˜g_{δ}(t_{2,i}+τ)−g˜_{δ}(t_{2,i}−τ)
.
and obtain a set of values {α_{i}, β_{i}} (i= 0, N);

• Step 4: Calculate the mean values α_{e} and β_{e}, which are taken as the solution of
the inverse problem.

In numerical experiments, t_{1,0} = τ_{1}, t_{2,0} = τ_{1,0} + 20τ. Further each successive pair
of (t_{1,i}, t_{2,i}) is obtained from the previous one by shifting one step τ = 1 (day) to the
right, N = 150.

Recovery results:

• Smoothing by moving average: recovered values – α_{e}= 3.18·10^{−2} andβ_{e} = 5.86;

the relative error – 6 and 17 % respectively;

• Smoothing by Gaussian kernel: recovered values –α_{e} = 3.18·10^{−2} andβ_{e} = 5.86;

the relative error – 6 and 17 % respectively.

Fig. 2 presents functionsg(t),g_{δ}(t)and ˜g_{δ}(t), as well as the recovery resultsα_{i} and
β_{i} (i= 0, N), with different smoothing options.

Smoothing with a Gaussian kernel gives a much smaller variation of values of α_{i}
and β_{i} relative to the exact values. However, in both cases, the average values α_{e} and
βe are identical.

The disadvantage of the proposed method is that its implementation assumes nu- merical differentiation of data which is an ill-posed problem [39].

### 4.2 Second numerical method of solving

To reduce the influence of errors contained in g˜δ(t), we modify the above-stated algo- rithm:

Figure 2: The functions g(t), g_{δ}(t), and g˜_{δ}(t); the results of recovery α_{i} and β_{i} (i =
0, N). a) and d) g(t) – solid line; g_{δ}(t) – dotted line; a) ˜g_{δ}(t) – solid thin line, data
smoothing by moving average; d)˜g_{δ}(t)– solid thin line, data smoothing by the Gaussian
kernel. b,c) the results of the recovery of values of α_{i} and β_{i} for different times t_{1,i}
and t_{2,i} (i = 0, N) in the interval [τ_{1}, τ_{2}], using data smoothing by moving average,
recovered values are shown by the dotted line, exact values – by the solid one. e,f) the
results of the recovery of values of α_{i} and β_{i} for different times t_{1,i} and t_{2,i} (i = 0, N)
in the interval [τ_{1}, τ_{2}], using data smoothing by the Gaussian kernel, recovered values
are shown by the dotted line, exact values – by the solid one.

Figure 3: Results of recovery of α_{i} and β_{i} (i = 0, N). Recovered values are shown by
the dotted line, the exact value – by the solid one.

• Step 1: Select the interval[τ_{1}, τ_{2}](see. Fig. 1), from which we take data to recover
the parameters α and β;

• Step 2: Choose time pointst_{1,i}< t_{2,i}andt_{3,i}< t_{4,i},t_{1,i}< t_{3,i},t_{2,i} < t_{4,i}, (i= 0, N)
from the interval [τ_{1}, τ_{2}];

• Step 3: Solve the system of equations

t2,i

Z

t1,i

˜

g_{δ}(t)dt·α−

t2,i

Z

t1,i

K(t)dt·β = ˜g_{δ}(t_{1,i})−g˜_{δ}(t_{2,i}),

t4,i (15) Z

t3,i

˜

g_{δ}(t)dt·α−

t4,i

Z

t3,i

K(t)dt·β = ˜g_{δ}(t_{3,i})−g˜_{δ}(t_{4,i}).

and get a set of values {α_{i}, β_{i}} (i= 0, N);

• Step 4: Calculate the mean values α_{e} and β_{e}, which are taken as the solution of
the inverse problem.

In system (15), the integrals are calculated using the trapezoidal rule. In numerical
experiments t_{1,0} = τ_{1}, t_{2,0} = t_{1,0} + 70τ, t_{3,0} = τ_{1}+ 50τ, t_{4,0} = t_{3,0} + 70τ. Every next
pair (t1,i, t2,i) and (t3,i, t4,i) are obtained from the previous ones by shifting one step
τ = 1 (day) to the right, N = 150.

Since the smoothing with a Gaussian kernel has given smaller spread of obtained
values, only this procedure was used for smoothing the functions g_{δ}(t).

Recovery results: αe = 3.07·10^{−2} and βe = 5.22, the relative error –2.36and 4.36%
respectively.

Fig. 3 shows the results of the recovery of valuesα_{i} and β_{i}. We see that an error of
the recovery has decreased.

It should be noted that the calculations for each pair of(t1,i, t2,i)and (t3,i, t4,i) can be conducted in parallel.

### 4.3 Third numerical method of solving

The following method can be considered as a traditional one. We write a function of the discrepancy:

J(α, β) =

N

X

i=1

n^{2}(t_{i}), n(t_{i}) = g(t_{i})−˜g_{δ}(t_{i}). (16)
To find the minimum we use Newton’s method (see, for example, [40]):

α^{k+1}
β^{k+1}

=
α^{k}

β^{k}

−H^{−1}(α^{k}, β^{k})∇J(α^{k}, β^{k}).

In order to use the method, we need to know the expression for the gradient and the matrix of second derivatives of the function J, which may be obtained as follows:

∇J(α, β) =
J_{α}

J_{β}

, H(α, β) =

J_{αα} J_{αβ}
J_{αβ} J_{ββ}

, where

J_{α} =

N

X

i=1

2n(t_{i})g_{α}(t_{i}), J_{β} =

N

X

i=1

2n(t_{i})g_{β}(t_{i}),

g_{α}(t) = −

t

Z

0

g(s)e^{−α(t−s)}ds, g_{β}(t) = −

t

Z

0

K(s)e^{−α(t−s)}ds,

J_{αα} =

N

X

i=1

2(g^{2}_{α}(t_{i}) +n(t_{i})g_{αα}(t_{i})), J_{ββ} =

N

X

i=1

2g_{α}^{2}(t_{i}),

J_{αβ} =

N

X

i=1

2(g_{α}(t_{i})g_{β}(t_{i}) +n(t_{i})g_{αβ}(t_{i})),

g_{αα}(t) = −2

t

Z

0

g_{α}(s)e^{−α(t−s)}ds, g_{αβ}(t) = −

t

Z

0

g_{β}(s)e^{−α(t−s)}ds.

It is easy to see that the functions introduced above are the solutions of the following problems:

g^{0}_{α}(t) +αg_{α}(t) = −g(t),
g_{α}(0) = 0.

g^{0}_{β}(t) +αg_{β}(t) = K(t),
g_{β}(0) = 0.

g^{0}_{αα}(t) +αg_{αα}(t) = −2g_{α}(t),
gαα(0) = 0.

g^{0}_{αβ}(t) +αg_{αβ}(t) = −g_{β}(t),
g_{αβ}(0) = 0.

The dependence of J of the arguments α and β is nonlinear. In general, this function may have local minima and maxima, thus the minimization process can lead to a local minimum. It is also known that Newton’s method requires that the initial

Figure 4: The shape of the function J(α, β). It reaches minimum at (α, β) = (3.0·
10^{−2},5). a) separately highlighted region J(α, β) < 0.1, b) separately highlighted
region J(α, β)<0.01

approximation is close enough to the minimum point. However, this problem could be resolved as follows (see., for example [40]): first, several steps are made using, for example, the Conjugate Gradient Method, and then, the Newton’s method is applied.

Since the approach is well known, we will not mention it onwards.

Fig. 4 shows the function J(α, β) which has a convex and ravine graph. Fig. 4a displays a mapped area where J[α, β]≤0.1, while Fig. 4b –J[α, β]≤0.01.

The figure shows that, on the one hand, the function has no local minima and maxima; on the other hand, it has a ravin (level lines extend in the same direction) and a sloping form (a region where the value of J is less than some ε, is sufficiently large). If there are errors in the data of the inverse problem then a region of the function values, where it is less than some ε, increases (see. [41, 42, 43]), which causes an error in searching for a global minimum. , It follows from the foregoing that the result of the minimization of the function J(α, β) may depend on the location of an initial approximation.

In order to apply the Newton’s method for the numerical solution of the inverse problem, we choose four different initial approximations. After minimizftion of function (16) we get the following recovery results:

• initial approximation – α^{0} = 1.0·10^{−2} and β^{0} = 1.0; result of minimization –
α^{∗} = 2.91·10^{−2} and β^{∗} = 4.72; relative error –3.0 and 5.6% respectively;

• initial approximation – α^{0} = 5.0·10^{−2} and β^{0} = 9.0; result of minimization –
α^{∗} = 3.34·10^{−2} and β^{∗} = 5.48; relative error –6.8 and 9.6% respectively;

• initial approximation – α^{0} = 1.0·10^{−2} and β^{0} = 9.0; result of minimization –
α^{∗} = 3.07·10^{−2} and β^{∗} = 5.18; relative error –2.33 and 3.6 % respectively;

• initial approximation – α^{0} = 5.0·10^{−2} and β^{0} = 1.0; result of minimization –
α^{∗} = 3.08·10^{−2} and β^{∗} = 4.81; relative error –2.67 and 3.8 % respectively.

### 4.4 Conclusions from numerical experiments

After trying three methods for solving inverse problem (1)-(5) we found that the first method is the least resistant to the presence of errors in the data. The reason for this is that the approach involves numerical differentiation, which is an ill-posed problem.

The second and third methods with an appropriate choice of the initial approximation give results comparable in accuracy of recovery. However, in terms of ease of imple- mentation, the speed of computation and possibilities for parallelization, the second approach to the numerical solution of the inverse problem is preferred.

### 4.5 Numerical solution of the inverse problem, the case of several ’humps’

All the above calculations were performed when data of the inverse problem had one

’hump’, i.e. the contamination has taken place, and it was possible to allocate a time
period when the given concentration C(x_{0}, t) has taken the form as, for example, in
Fig. 1. However, there could be a case when at x = 0, at some time intervals several
facts of groundwater contamination occur. In this case, the curve g(t) may have two
or more ’humps’, see Fig. 5.

Figure 5: An example of long-term observations of the concentration of a contaminant in groundwater. a)-c) – stages of the recovery of parameters associated with each

’hump’.

• after solving direct problem (12) with the help of found α and β, obtain the function g(t) and subtract it from the long-term records data (see. Fig. 5a);

• define a new initial time τ_{0} (see Fig. 5d), which will be considered as the new
starting point of countdown;

• select a new time interval [τ_{1}, τ_{2}] which will be used for the recovery of the pa-
rameters α_{2} and β_{2} for the second ’hump’;

• after solving direct problem (12) with the help of found α and β, obtain the function g(t) and subtract it from the long-term records data (see Fig. 5c);

• re-define the initial time τ_{0} (see Fig. 5с), which will be taken as the new starting
countdown;

• select a new time interval [τ_{1}, τ_{2}] which will be used for the recovery of the pa-
rameters α_{3} and β_{3} for the third ’hump’;

For solving of inverse problems of determining the parameters α_{k} and β_{k} we used
the second method. Recovery results are below:

• exact values – α_{1} = 3.0·10^{−2} and β_{1} = 5.0; recovered values – α_{1e} = 2.93·10^{−2}
and β_{1e}= 5.28; relative error – 2.33 and 5.6% respectively;

• exact values – α_{2} = 2.0·10^{−2} and β_{2} = 4.0; recovered values – α_{2e} = 2.09·10^{−2}
and β2e= 4.28; relative error – 4.5 and 7.0% respectively;

• exact values – α_{3} = 5.0·10^{−2} and β_{3} = 7.0; recovered values – α_{3e} = 5.15·10^{−2}
and β_{3e}= 7.31; relative error – 3.0 and 4.43% respectively.

Thus, we have quite satisfactory recovery results. However, in the example above, you can clearly see the presence of three ’humps’, i.e. there were three facts of ground- water contamination spaced over time. The second contamination took place 120 days after the first, and the third – a year later.

Fig. 6a displays an example of the function g(t), when the second contamination took place 10 days after the first, and the third – 25 days. We can see that all the

’humps’ have merged. We apply the recovery method to the data in Fig. 6 assuming that only one fact of pollution took place, i.e. the contamination is under conditions of model (11). It is reasonable to assume that, in this case, the recovered parameter α would be close to the lowest of three αk, while β – to the sum of the three βk, but slightly less since the facts of contamination are separated in time.

The result of recovery isα_{e}= 2.48·10^{−2} and β_{e}= 14.11. Fig. 6a and b presents the
results of calculation of α_{i} andβ_{i}, where data of the inverse problem were taken either
without noise or with noise P = 3%. The fact that the ’hump’ shown in Fig. 6a is
the sum of several ’humps’ is confirmed by the presence of a curved trend component
in Fig. 6a and b (compare with Fig. 2b,c,e and f, and Fig. 3). The presence of the
trend component tells us that the adopted model of contamination (11) does not work
properly and should be replaced by another one.

If we apply the third method for the recovery of the parameters α and β (the
method of minimizing of the target functional), we get the curve g(t)that is very close
to the curve in Fig. 6a and some values α^{∗} and β^{∗}, which will be taken as a solution
of the inverse problem (in the numerical simulations, we obtained: α^{∗} = 2.38·10^{−2} and
β^{∗} = 14.03). However, we will not have any information that the adopted model of
contamination (11) is not suitable for the interpretation.

### 5 The recovery in the case of merging several ’humps’

As we have already noted, the theorem of existence of solutions of the Volterra integral
equation of the first kind (10), which can not be reduced to the Volterra integral
equation of the second kind, in general, is not proved. However, in a particular case, it
is possible to formulate conditions for the existence of a solution of integral equation
(10) in the form of (11). The proof was constructive suggesting a method of finding
of a solution of integral equation (10). Nevertheless, the results of the numerical
calculations have shown that the model of a single source of contamination (11) could
be inapplicable because there could be several sources of type (11) little spaced over
time. The responses from them at a distance x=x_{0} could merge and the result would
be similar to the response from a single pollution source. The presence of the trend

Figure 6: a) An example of observations of the concentration of a contaminant in groundwater, when the ’humps’ of three different facts of contamination have merged;

b,с) the results of the recovery values α_{i} and β_{i}, the solid line – the recovery is based
on the accurate data, the dotted line – the recovery is based on data with noise.

assuming that the right part of g(t)satisfies the condition g(0) = 0.

We do not assume that any of the conditions allowing to reduce the integral equation
(10) to the Volterra integral equation of the second kind are met. On the contrary, for
example, a case is possible when K^{(n)}(0) = 0 (n = 0,1,2...) that does not allow such
reducing.

Consider a spaceF^{β}[0, T](β >1) of continuous functionsf(t)on the interval[0, T]
decomposable in a Fourier series

f(t) = ˆb_{0}+

∞

X

k=1

ˆ

a_{k}sin(ω_{k}t) + ˆb_{k}cos(ω_{k}t)

, ω_{k} = 2πk/T, (17)

and |ˆa_{k}| ≤Ak^{−β} and |ˆb_{k}| ≤Bk^{−β} (b_{0}, A >0and B >0 are some constants).

Denote O_{0}^{ε} to be a punctured ε-neighborhood of the point t = 0.

Theorem 1. Suppose that K(t) 6≡ 0 is a piecewise continuous function on the
interval [0, T] and |K(t)| ≤ m_{0}. A solution f(t)∈F^{β}[0, T] (β = 1 +α, α > 0) of the
integral equation (1) exists if and only if, on the interval [0, T], the function g(t), can
be represented as

g(t) =b_{0}g_{0}^{0}(t) +

∞

X

k=1

a_{k}g_{k}(t) +b_{k}g^{0}_{k}(t)

, (18)

where

1. functions gk(t) (k = 0,1, ...) are the solutions of problems
g_{k}^{00}(t) +ω_{k}^{2}g_{k}(t) = K(t),

g_{k}(0) = 0, g_{k}^{0}(0) = 0; (19)
2. |a_{k}| ≤Ak^{−α} and |b_{k}| ≤Bk^{−(1+α)}, and A >0, B >0 and b_{0} are some constants.

The solution takes the form

f(t) =b_{0}+

∞

X

k=1

a_{k}

ω_{k}sin(ω_{k}t) +b_{k}cos(ω_{k}t)

. (20)

Necessity. Let the solution of the integral equation (10) f(t)∈F^{1+α}[0, T], where
α > 0, then due to the uniform convergence of the series (17) the function g(t) is

representable in the form

g(t) = ˜b_{0}g_{0}^{0}(t) +

∞

X

k=1

˜a_{k}

t

Z

0

K(s) sin(ω_{k}(t−s))ds+ ˜b_{k}

t

Z

0

K(s) cos(ω_{k}(t−s))ds

.

It is obvious that the functions
g_{0}(t) =

t

Z

0

(t−s)K(s)ds, g_{k}(t) = 1
ω_{k}

t

Z

0

K(s) sin(ω_{k}(t−s))ds (21)

are solutions of the problem (19), and g_{k}(t)∈ C^{1}[0, T]. Thus, the function g(t) could
be written in the form (18) and b0 = ˜b0, ak = ˜akωk and bk = ˜bk (k = 1,2, ...). Since

|K(t)| ≤ m_{0}, the series converges uniformly. Hence, it converges to a continuous
function. Due to the boundary conditions (19) g(0) = 0.

Sufficiency. Let the function g(t)can be represented in the form (18). A solution of problem (19) could be represented in the form (21). It is easy to see that the following equality takes place:

t

Z

0

K(s)

"

b_{0}+

∞

X

k=1

a_{k}

ω_{k}sin(ω_{k}(t−s)) +b_{k}cos(ω_{k}(t−s))

#

ds=g(t).

We have |˜a_{k}|=|a_{k}/ω_{k}| ≤2πT^{−1}Ak^{−(1+α)} and|˜b_{k}|=|b_{k}| ≤Bk^{−(1+α)}, then the series in
brackets converges uniformly, therefore it converges to a continuous function f(t), and
f(t)∈F^{1+α}[0, T].

If we assume that there exists such ε > 0 that K(t) 6= 0 when t ∈ O_{0}^{ε}, then
the uniqueness of the solution of integral equation (10) follows from the Titchmarsh
Theorem [35] (see also [36, стр. 49]).

Theorem 2. Let there exist two functions, g^{1}(t) and g^{2}(t), for which following
representations are valid

g^{i}(t) =b^{i}_{0}g_{0}^{0}(t) +

∞

X

k=1

a^{i}_{k}g_{k}(t) +b^{i}_{k}g^{0}_{k}(t)

, i= 1,2,

where b^{i}_{0}, a^{i}_{k}andb^{i}_{k}(k = 1, k,i= 1,2) satisfy condition 2 from Theorem 1. IfK(t)6≡0,

|b^{1}_{0} −b^{2}_{0}| ≤ δ, |a^{1}_{k}−a^{2}_{k}| ≤ δk^{−α}, |b^{1}_{k}−b^{2}_{k}| ≤ δk^{−(1+α)} (k = 1,2, ...,) then the solutions
f^{1}(t) and f^{2}(t) of integral equation (1) satisfy the following inequality:

|f^{1}(t)−f^{2}(t)| ≤Cδ,
where C =C(T, α).

First of all, it is obvious that the difference g(t) = g^{1}(t)−g^{2}(t) could be repre-
sented in form (18). Hence, the solution f(t) = f^{1}(t)−f^{2}(t) of equation (1), which

= 1 + (T /2π+ 1)ζ(1 +α) δ≡C(T, α)δ, where ζ(·) is Riemann function.

Note that Theorem 2 implies the uniqueness and continuous dependence on the
input data of the solution of integral equation (10) inF^{β}[0, T]without the assumption
that K(t)6= 0 att ∈ O^{ε}_{0}.

### 5.2 Reconstruction method

Based on the results of the previous section, we search the function f(t)approximating the function g(t)by the sum

b0g_{0}^{0}(t) +

N

X

k=1

akgk(t) +bkg_{k}^{0}(t)

by means of minimization of the residual functional

J[b_{0}, ..., a_{k}, ..., b_{k}, ...] =

N

X

j=0

"

g(t_{j})−b_{0}g_{0}^{0}(t_{j})−

N

X

k=1

a_{k}g_{k}(t_{j})+b_{k}g_{k}^{0}(t_{j})

#^{2}

. (22) The number N in both sums is selected based on the following considerations. If, for example, observations are being made during one year, while the concentration of pollution data are taken once a day, then N = 365.

Residual functional (22) is quadratic and for its minmization the Conjugate Gradi- ent Method is applicable (see, for example, [40]).

We consider the example from Section 4.4 with ’humps’ merged, i.e. the function g(t) has the form shown in Fig. 6a. We minimize residual functional (22) using the Method of Conjugate Gradients (due to its evidence, we do not provide the expression for the gradient of the residual functional). The results of calculations are presented in Fig. 7 and data error g(t)is 3%.

As already mentioned, we model the pollution by a decaying exponential because its coefficients and index have physical meaning. A practician needs them for the interpretation of the results. Looking at the form of the restored curve g(t) in Fig. 7 we search the function f(t) as follows

f(t) =

3

X

m=1

β_{m}e^{−α}^{m}^{(t−t}^{m}^{)}θ(t−t_{m}), (23)

Figure 7: The results of restoring of the function f(t). The exact value is highlighted by the solid line, the restored one - by the dotted line.

where θ(t) is Heaviside function. In this case, the number of contamination cases is
equal to three and we can roughly estimate valuesα_{m}, β_{m},t_{m},t = 1,3,t_{1} = 0. For the
more accurate values we minimize the residual functional:

J[..., α_{m}, ..., β_{m}, ..., t_{m}, ...] =

N

X

j=0

g(t_{j})−

tj

Z

0

K(t_{j} −s)f(s)ds

2

. (24)

Recovery results:

• exact values – α_{1} = 3.0·10^{−2} and β_{1} = 5.0; recovered values – α_{1e} = 2.94·10^{−2}
and β_{1e}= 5.11;

• exact values – α2 = 2.0 · 10^{−2}, β2 = 4.0, and t2 = 10; recovered values –
α_{2e}= 2.10·10^{−2}, β_{2e}= 4.23, and t_{2} = 10;

• exact values – α_{3} = 5.0 · 10^{−2}, β_{3} = 7.0, and t_{2} = 25; recovered values –
α_{3e}= 5.17·10^{−2}, β_{3e}= 7.29, and t_{2} = 26.

It should be noted that if we skip the first step - the minimization of the residual functional (22) - and immediately seek a solution as sum (23), it is not clear how many terms of the sum to take. In addition, after the minmimization of the residual functional (22), approximate values of the unknown parameters can be obtained. It is also worth mentioning that the proposed method of finding a solution of integral equation (10) used at the first step of determining the parameters of the source of groundwater pollution could be applied in general case of the solution of the Volterra integral equation of the 1st kind of the convolution type. To the best of our knowledge, this method is not mentioned in the reference literature [33, 44, 45].

computing speed and a potential for parallelization.

Where the facts of contamination were sufficiently separated in time from each other and appeared in the observed function g(t)in the form of a separate ’hump’ the proposed algorithm allowed to determine the parameters of the source of contamination rather well. However, when the facts of contamination were separated in time insuf- ficiently and merged into one ’hump’ of function g(t), the proposed method produced some integral characteristics of the pollution parameters and an indication that the chosen model of contamination is not suitable. In the case where the ’humps’ of the function g(t) were merged into one, a method based on the minimization of a residual functional have been also proposed allowing to calculate the parameters of a pollution source in this general case.

### Acknowledgement

The work of the 1th author was supported by the project of the program of Presidium of RAS № 0314-2015-0011. The work of the 2nd, 3nd, and 4th authors was supported by the Project “Development of information-analytical system of monitoring of the reserves and quality of groundwater of the Republic of Kazakhstan” (80-420-16).

### References

[1] Rogachevskaya L.M., Estimation of water exchange rate of groundwater north-eastern part of the Dnieper artesian basin, Water resources, Vol. 27, n. 5 (2000), p. 566-573. (In Russian) [2] Rogachevskaya L.M., The study of migration features of Chernobyl of Cs-137 in the aeration

zone for identifying regional vulnerabilities to groundwater radionuclide contamination, In the book: Modern problems of hydrogeology and hydromechanics, St. Petersburg, 2002, p. 146-150.

(In Russian)

[3] Turganbaev E., Nugumanova A., Design of the remote access to the “Ground water of the Re- public of Kazakhstan” database, Proceedings of International Scientific and Practical Conference

“Green economy is the future of humanity”, 2014, p. 1665-1677.

[4] Foster S.,Fundamental conception aquifer vulnerability, pollution risk and ground water protec- tion, TNO Committee of Hydrogeological Research, The Hange, Proceeding and Information, n.

38 (1987), p. 69-86.

[5] Schreiber R.P., The National Ground Water Monitoring Network Data Portal: From Pilot to Production, NGWA Summit, Proceedings of The National and International Conference on Groundwater, NGWA, 2013, p. 212-219.

[6] Van Dijk A., Renzullo L.J., Water resource monitoring systems and the role of satellite obser- vations, Hydrology and Earth System Sciences, Vol. 15, n. 1 (2011), p. 39-55.

[7] Panichkin V.Y., Miroshnichenko O.L., Trushel L.Y., Zakharova N.M., Automated information system “Groundwater of Kazakhstan”.

[8] Beldeubayeva Zh., Rakhmetullina S., Turganbayev E., Saparbekov Zh., Effective Management of Groundwater System Based on Information Technology, Global Journal on Technology, n. 8 (2015), p. 245-253.

[9] Beldeubaeva Z.T., Rahmetullina S.Z., Turganbaev E.M., The concept of development of information-analytical system of ecological monitoring of groundwater, Vestnik of EKSTU, Ust- Kamenogorsk, n. 3 (2014), p. 145-149. (In Russian)

[10] Beldeubayeva Zh., Rakhmetullina S., Turganbayev E., Krivykh V., Information system of groundwater monitoring, Proceedings of 15th International multidisciplinary scientific geocon- ference SGEM 2015, Bulgaria, Vol. 1 (2015), p. 139-146.

[11] Beldeubayeva Zh., Rakhmetullina S., Turganbayev E., Krivykh V.Information system of efficient data management of groundwater monitoring the Republic of Kazakhstan, Proceedings of 9th International Conference on Application of Information and Communication Technologies AICT 2015, Russia, 2015, p. 195-201.

[12] Mikhailov V.N., Dobrovolsky A.D., Dob S.A.,Hydrology: The textbook for universities, 3rd ed., Moskva, Vysshaya shkola, 2008.

[13] Bochever F.M., Lapshin N.N., Oradovskaya A.E.,Protection of groundwater against pollution, Moskva, Nedra, 1979.

[14] Javandel I., Doughty I., Tsang C.-F.,Groundwater transport: Handbook of mathematical models, American Geophysical Union, Water resources monograph series, vol. 10, Washington, D.C., 1984.

[15] Atmadja J., Bagtzoglou A.C.,State of the Art Report on Mathematical Methods for Groundwater Pollution Source Identification, Environmental Forensics, n. 2 (2001), p. 205-214.

[16] Bhattacharjya R.K., Solving Groundwater Flow Inverse Problem Using Spreadsheet Solver, J.

Hydrol. Eng., Vol. 16, n. 5 (2011), p. 472-477.

[17] Carrera J., State of the Art of the Inverse Problem Applied to the Flow and Solute Transport Equations, In the book: Groundwater Flow and Quality Modelling, NATO ASI Series, Springer Netherlands, 1988, p. 549-583.

[18] Carrera J., Medina A., Galarza G.,Groundwater inverse problem. Discussion on geostatistical formulations and validation, Hydrogeologie, n. 4 (1993), p. 313-324.

[19] Ginn, T.R., Cushman, J.H.,Inverse methods for subsurface flow: A critical review of stochastic techniques, Stochastic Hydrology and Hydraulics, Vol. 4 (1990), p. 1-26.

[20] Huang Ch.-H., Li J.-X., Kim S.,An inverse problem in estimating the strength of contaminant source for groundwater systems. Applied Mathematical Modelling, Vol. 32 (2008), p. 417–431.

[21] Kitanidis P, Vomvoris E.G.,A geostatistical approach to the problem of groundwater modelling (steady state) and one-dimensional simulation, Water Resources Research, 1983, v. 19, n. 3, p.

677-690.

[22] Kriksin Y.A., Ivy S.N., Samarskaya E.A., Tishkin V.F., The inverse problem of source recon- struction for convective diffusion equation. Matem. modelirovanie, Vol. 7, n. 11 (1995), p. 95-108.

resources reseach, Vol. 32, n. 5 (1996), p. 1131–1161.

[26] Nguen V.T., Nguen H.T., Tran T.B., Vo A.K.,On an inverse problem in the parabolic equation arising from groundwater pollution problem, Boundary Value Problems, Vol. 67 (2015), p. 1-23.

[27] Sun N.-Zh., Inverse Problems in Groundwater Modeling Series: Theory and Applications of Transport in Porous Media, Vol. 6, Kluwer Academic Publishers, Dordrecht, 1994.

[28] Sun N.-Zh., Yeh W.W.-G., Coupled Inverse Problems in Groundwater Modeling 1. Sensitivity Analysis and Parameter Identification, Water Resources Research, Vol. 26, n. 10 (1990), p.

2507-2525.

[29] Woodbury A.D., A full-Bayesian approach to the groundwater inverse problem for steady state flow, Water resources research, Vol. 36, n. 8 (2000), p. 2081–2093.

[30] Yeh W.W.-G., Review of parameter identification procedures in groundwater hydrology: The inverse problem, Water Resources Research, Vol. 22, n. 2 (1986), p. 95-108.

[31] Alifanov O.M., Artyukhin E.A., Rumyantsev S.V.,Extreme Methods for Solving Ill-Posed Prob- lems Extreme Methods for Solving Ill-Posed Problems with Application to Inverse Heat Transfer Problems, Begell House Inc., New York, 1996.

[32] Bateman H., Erdelyi A., Tables of Integral Transforms, Vol. 1, McGraw-Hill Book Co., New York, 1954.

[33] Manzhirov A.V., Polyanin A.V.,Handbook of Integral Equations: methods of solution, Moskva, Factorial Press, 2000.

[34] Tikhonov A.N., Arsenin V.J., Methods of solving ill-posed problems, Moskva, Nauka, 1979.

[35] Titchmarsh E.C.,The Theory of Functions, Oxford University Press, 1939.

[36] Bukhgeim A.L.,Volterra equations and inverse problems, Novosibirsk, Nauka, 1983.

[37] Doudkin V.A., Akimova Y.S.,Research of methods of smoothing the seismic signal, Proceedings of the International Symposium “Safety and quality”, n. 1 (2005), p. 367-370.

[38] Numerical Recipes in Fortran 77: The Art of Scientific Computing, Cambridge University Press, 1992.

[39] Lavrentiev M.M., Romanov V.G., Shishatskiy S.P.,Incorrect Problems of mathematical physics and analysis, Moskva, Nauka, 1980.

[40] Vasilyev F.P.,Numerical methods for solving extreme problems, Moskva, Nauka, 1988.

[41] Karchevsky A.L.,Simultaneous reconstruction of permittivity and conductivity, Journal of Inverse and Ill-Posed Problems, Vol. 17, n. 4 (2009), p. 385-402.

[42] Nazarova, L.A., Nazarov, L.A., Karchevsky, A.L., Vandamme, M.,Estimating diffusion-capacity parameters of a coal bed using the gas pressure measured in a hole and the solution of an inverse problem, Journal of Applied and Industrial Mathematics, Vol. 8, n. 2 (2014), p. 267-273.

[43] Duchkov A.A., Karchevsky A.L., Determination of Terrestrial Heat Flow from Temperature Measurements in Bottom Sediments, Journal of Applied and Industrial Mathematics, Vol. 7, n.

4 (2013), p. 480-502.

[44] Verlan’ A.F., Sizikov V.S., Integral equations: methods, algorithms, codes. Handbook, Kiev, Naukova dumka, 1986, 543 p. (In Russian)

[45] EqWorld, http://eqworld.ipmnet.ru

Karchevsky A.L.

Sobolev Institute of Mathematics SB RAS, Russia, 630090 Novosibirsk, pr. Koptyuga, 4 Email: karchevs@math.nsc.ru

Turganbayev Y.M.

Graduate School of Business, Nazarbayev University,

Kazakhstan, 010000 Astana, Kabanbay Batyr Ave, 53, block C3 Email: yerken.turganbayev@gmail.com,

Rakhmetullina S.J.

D. Serikbayev East Kazakhstan State Technical University, Kazakhstan, 070004 Ust-Kamenogorsk, st. A.K. Protozanova, 69 Email: SRakhmetullina@ektu.kz,

Beldeubayeva Zh.T.

D. Serikbayev East Kazakhstan State Technical University, Kazakhstan, 070004 Ust-Kamenogorsk, st. A.K. Protozanova, 69 Email: leka1997@mail.ru

Received 11.10.2016, Accepted 01.11.2016