Numerical solution of an inverse problem of determining the parameters of a source of groundwater pollution

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₀ 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₀, 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₀:

C(x₀, 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ˆ ₀, p) = ˆg(p). (8)

We find the solution of equation (6) in general form C(x, p) =ˆ A₁e^λ1x+A₂e^λ2x, where λ_1,2 = 1

2

"

ν D ±

rν²

D² + 4R D

# .

According to the second boundary condition (7) A₁ = 0. (8) leads to C(x, p) = ˆˆ g(p)e^λ2(x−x0).

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²+p

a+p a²+p

f(p),ˆ where a= ν

√RD, b=x₀ 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+a2terfc b

2√

t +a√ t

e^−a2t.

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³/d) characterize the intensity of groundwater contamination at x= 0.

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

g⁰(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⁻², β = 5.0.

Data error of the inverse problem was P = 3%. D = 10(m²/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^−t2(t), (13)

where

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

η(τ) = −1.26551223 + 1.00002368τ+ 0.37409196τ²+ 0.09678418τ³

−0.18628806τ⁴+ 0.27886807τ⁵−1.13520398τ⁶ + 1.48851587τ⁷

−0.82215223τ⁸+ 0.17087277τ⁹.

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 [τ₁, τ₂];

• 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 = τ₁, 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⁻² andβ_e = 5.86;

the relative error – 6 and 17 % respectively;

• Smoothing by Gaussian kernel: recovered values –α_e = 3.18·10⁻² 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 [τ₁, τ₂], 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 [τ₁, τ₂], 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[τ₁, τ₂](see. Fig. 1), from which we take data to recover the parameters α and β;

• Step 2: Choose time pointst_1,i< t_2,iandt_3,i< t_4,i,t_1,i< t_3,i,t_2,i < t_4,i, (i= 0, N) from the interval [τ₁, τ₂];

• 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 = τ₁, t_2,0 = t_1,0 + 70τ, t_3,0 = τ₁+ 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⁻² 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²(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⁻¹(α^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²_α(t_i) +n(t_i)g_αα(t_i)), J_ββ =

N

X

i=1

2g_α²(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⁰_α(t) +αg_α(t) = −g(t), g_α(0) = 0.

g⁰_β(t) +αg_β(t) = K(t), g_β(0) = 0.

g⁰_αα(t) +αg_αα(t) = −2g_α(t), gαα(0) = 0.

g⁰_αβ(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⁻²,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 – α⁰ = 1.0·10⁻² and β⁰ = 1.0; result of minimization – α^∗ = 2.91·10⁻² and β^∗ = 4.72; relative error –3.0 and 5.6% respectively;

• initial approximation – α⁰ = 5.0·10⁻² and β⁰ = 9.0; result of minimization – α^∗ = 3.34·10⁻² and β^∗ = 5.48; relative error –6.8 and 9.6% respectively;

• initial approximation – α⁰ = 1.0·10⁻² and β⁰ = 9.0; result of minimization – α^∗ = 3.07·10⁻² and β^∗ = 5.18; relative error –2.33 and 3.6 % respectively;

• initial approximation – α⁰ = 5.0·10⁻² and β⁰ = 1.0; result of minimization – α^∗ = 3.08·10⁻² 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₀, 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 τ₀ (see Fig. 5d), which will be considered as the new starting point of countdown;

• select a new time interval [τ₁, τ₂] which will be used for the recovery of the pa- rameters α₂ and β₂ 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 τ₀ (see Fig. 5с), which will be taken as the new starting countdown;

• select a new time interval [τ₁, τ₂] which will be used for the recovery of the pa- rameters α₃ and β₃ 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 – α₁ = 3.0·10⁻² and β₁ = 5.0; recovered values – α_1e = 2.93·10⁻² and β_1e= 5.28; relative error – 2.33 and 5.6% respectively;

• exact values – α₂ = 2.0·10⁻² and β₂ = 4.0; recovered values – α_2e = 2.09·10⁻² and β2e= 4.28; relative error – 4.5 and 7.0% respectively;

• exact values – α₃ = 5.0·10⁻² and β₃ = 7.0; recovered values – α_3e = 5.15·10⁻² 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⁻² 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⁻² 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₀ 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⁽ⁿ⁾(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₀+

∞

X

k=1

ˆ

a_ksin(ω_kt) + ˆb_kcos(ω_kt)

, ω_k = 2πk/T, (17)

and |ˆa_k| ≤Ak^−β and |ˆb_k| ≤Bk^−β (b₀, A >0and B >0 are some constants).

Denote O₀^ε 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₀. 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₀g₀⁰(t) +

∞

X

k=1

a_kg_k(t) +b_kg⁰_k(t)

, (18)

where

1. functions gk(t) (k = 0,1, ...) are the solutions of problems g_k⁰⁰(t) +ω_k²g_k(t) = K(t),

g_k(0) = 0, g_k⁰(0) = 0; (19) 2. |a_k| ≤Ak^−α and |b_k| ≤Bk^−(1+α), and A >0, B >0 and b₀ are some constants.

The solution takes the form

f(t) =b₀+

∞

X

k=1

a_k

ω_ksin(ω_kt) +b_kcos(ω_kt)

. (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₀g₀⁰(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₀(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¹[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₀, 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₀+

∞

X

k=1

a_k

ω_ksin(ω_k(t−s)) +b_kcos(ω_k(t−s))

#

ds=g(t).

We have |˜a_k|=|a_k/ω_k| ≤2πT⁻¹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₀^ε, 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¹(t) and g²(t), for which following representations are valid

gⁱ(t) =bⁱ₀g₀⁰(t) +

∞

X

k=1

aⁱ_kg_k(t) +bⁱ_kg⁰_k(t)

, i= 1,2,

where bⁱ₀, aⁱ_kandbⁱ_k(k = 1, k,i= 1,2) satisfy condition 2 from Theorem 1. IfK(t)6≡0,

|b¹₀ −b²₀| ≤ δ, |a¹_k−a²_k| ≤ δk^−α, |b¹_k−b²_k| ≤ δk^−(1+α) (k = 1,2, ...,) then the solutions f¹(t) and f²(t) of integral equation (1) satisfy the following inequality:

|f¹(t)−f²(t)| ≤Cδ, where C =C(T, α).

First of all, it is obvious that the difference g(t) = g¹(t)−g²(t) could be repre- sented in form (18). Hence, the solution f(t) = f¹(t)−f²(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^ε₀.

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₀⁰(t) +

N

X

k=1

akgk(t) +bkg_k⁰(t)

by means of minimization of the residual functional

J[b₀, ..., a_k, ..., b_k, ...] =

N

X

j=0

"

g(t_j)−b₀g₀⁰(t_j)−

N

X

k=1

a_kg_k(t_j)+b_kg_k⁰(t_j)

#²

. (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

β_me^{−αm(t−tm)}θ(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₁ = 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 – α₁ = 3.0·10⁻² and β₁ = 5.0; recovered values – α_1e = 2.94·10⁻² and β_1e= 5.11;

• exact values – α2 = 2.0 · 10⁻², β2 = 4.0, and t2 = 10; recovered values – α_2e= 2.10·10⁻², β_2e= 4.23, and t₂ = 10;

• exact values – α₃ = 5.0 · 10⁻², β₃ = 7.0, and t₂ = 25; recovered values – α_3e= 5.17·10⁻², β_3e= 7.29, and t₂ = 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: [email protected]

Turganbayev Y.M.

Graduate School of Business, Nazarbayev University,

Kazakhstan, 010000 Astana, Kabanbay Batyr Ave, 53, block C3 Email: [email protected],

Rakhmetullina S.J.

D. Serikbayev East Kazakhstan State Technical University, Kazakhstan, 070004 Ust-Kamenogorsk, st. A.K. Protozanova, 69 Email: [email protected],

Beldeubayeva Zh.T.

D. Serikbayev East Kazakhstan State Technical University, Kazakhstan, 070004 Ust-Kamenogorsk, st. A.K. Protozanova, 69 Email: [email protected]

Received 11.10.2016, Accepted 01.11.2016