10th Asian Regional Conference of IAEG (2015)

Temperature-and-filtration conditions in the embankment at the Kumtor Mine,
Kyrgyzstan: Case study and contingency situation analysis
Larisa NAZAROVA(1), Leonid NAZAROV(1),Murataly JAMANBAEV(2) and Mirlan CHYNYBAEV(2)(1) Institute of Mining Siberian Branch of RAS, Novosibirsk, Russia
E-mail: lanazarova@ngs.ru(2) Razzakov Kyrgyz State Technical University, Bishkek, Kyrgyzstan
Abstract
The geomechanical model describing the heat-and-mass transfer processes in thevicinity of the protective embankment at the tailing pond of the Kumtor gold deposit(Kyrgyzstan) was designed based on the information on the structure of the object, as well as
the data on thermophysical properties and permeability of frozen and thawed ground andseasonal variation of air temperature. The simulation implemented with the help of original
FEM code has shown that: a) the zero isotherm between the frozen and thawed rocks attains
stable location 10–15 years after the embankment creation and the pond filling; b) the
decrease in the mean annual temperature essentially reduces infiltrate volume due to longer
time when impermeable frozen layer covers the embankment slope surface. The inverse
problem on determination of location and onset time of impervious screen damage by thedata of piezometric measurements was defined. Numerical experiments revealed that theformulated problem is unambiguously resolvable at an acceptable accuracy for any realheat-transfer and filtration properties of soil if input data were recorded at a number of points
and given the measurement error is moderate.
Keywords: protective embankment, heat, filtration, inverse problem
1.

Introduction

Being the largest in Central Asia, Kumtor gold
deposit was discovered in 1978 and put into operationin the mid-1990s. The deposit lies in the south-east
Kyrgyzstan, in the Tian Shan mountains, at anelevation of 4 km, in the permafrost zone
(www.kumtor.kg). This area belongs in the continental
climate zone RF Construction norms (1996), thewinter and summer temperatures are, respectively,– 35°C and + 20°C; the mean annual temperature is
6°C (Nepomnyashchaya, 2007).
Kumtor gold ore mine uses open pit mining
method, by now the mine has produced 5.6 Mt of ore
with the average gold content of 3.85 g/t(www.gold-deposit.ru/2183.html). Gold is recoveredby cyanidation, liquid waste goes to a tailing pond. As
the tailing pond fills up, its protective embankment isbuilt up by 1.5 m per year. At the present time, theembankment is 36 m high, the dam site is 280 m wide
and its extension Y=3 km (Osmonbetova, 2011). To
prevent slope erosion and infiltration, an impervious

Construction and operation of various purpose
waterworks (protection embankment, dams, dikes,
hydroelectric power plants) in permafrost zones
causes natural heat disbalance, rises temperature offrozen rocks and generates new migration paths for
fluid flow (Yershov and Williams, 2004; Chzhan,
2002; Vozin, 1987). This results in the instability of awaterworks structure and in pollution of environment
with production waste. Various aspects of modeling ofheat-mass exchange in terms of objects in permafrost
areas are discussed, for instance, in (Riseborough et al,2008; Yi et al, 2006). Substantial research deals withanalyzing behavior of large natural objects overgeological period of time (Bense and Person, 2008;
Rawlins et al., 2009; Velicogna et al., 2012; Zhang etal., 2008) and with studying behavior of production
and man-made objects (Goy et al., 1996; Lolaev et al.,
1997; Buiskikh and Zamoshch, 2010).

1

10th Asian Regional Conference of IAEG (2015)

Fig. 1 Profile of the object under investigation
norms, 1996). Let the embankment temperature andthe underlying ground temperature T = Ti at the initialtime moment t = 0, ant later on the tailing pond is filledwith the fluid with the temperature Tw , which isassumed constant in view of the continuity of the
production process.
Thawing is much slower than the heat exchange
between the fluid and the ground, therefore evolutionof the temperature T is described by the transient
equation of thermal conduction, considering
convection (Bejan, 2013):

screen has been created: the embankment is covered
with an impermeable 1.5 mm thick film that is sealed
along edging (www.vb.kg/244692). The weld seals are
the weakest place of the protection cover.
The temperature and filtration processes in the areaof the Kumtor mine embankment and the diagnostics
method for the state of the impervious screen are thesubjects of the present research.
2.

Formulation and solution of direct
problem

Figure 1 shows structural configuration of profileof the object under analysis. The tailing pond is
non-freezing over the entire service life, the fluid
temperature Tw varies from 4 to 8°C. The state of theprotective embankment and underlying ground ismonitored by means of measuring pressure andtemperature at different depth in observation
boreholes.
As the extension of the object is much more thanits height and width and shapes of the profiles of theembankment show a minor difference, the analysis
involves a two-dimensional problem in the Cartesian
coordinates (x, z). The current state of theembankment is of the prime importance, so we
simplify the situation by disregarding the change in theconfiguration of the embankment in the course ofoperation.
The lower boundary z = Z of the computational
domain (Fig. 2) lies in the zone of permafrost at thepermanent negative temperature Ti that is accepted toequal the temperature of –2°C of the neutral layer inthe discussed geographical region (RF Construction

T, t  V  grad T  div (kC grad T ) ,(1)

where kC = kC(T) is the thermal diffusivity coefficient;V is the fluid flow velocity defined by the Darcy law:V   k S grad P / pa ,(2)

kS = kS(T) is the filtration coefficient; pa is atmosphericpressure. Filtration takes place only in the thawing
zone (T > 0) where the additional pressure P obeys theequation:P, t  div (kS grad P) .

(3)

The boundary conditions are formulated asfollows (refer to Fig. 2):
at the section A0A1
T = Tw ,2

P = ρg(z(x) - z0)(4)

10th Asian Regional Conference of IAEG (2015)

Fig. 2 Computational domain, boundary conditions and finite element mesh fragment(ρ is the fluid density; g is the free fall acceleration;z(x) is the ordinate of current boundary point);
at A1…A8 (impervious screen) 
T=Tw , V  n  0

Table 1 Heat–filtration characteristics of rocks
Rock
106kS , m/s 106kC , m2/s
Thawed2.0
1.2
Frozen0
2.0(5)

Figure 3 demonstrates location of the zero
isotherm (the interface of the thawed and frozen
zones) at different time after filling the tailing pond at
Tw = 4°C (blue lines) and Tw = 8°C (red lines).

( n is the external normal to the boundary);at the daylight surface A8…A11
T=Ta , P=0,(6)

where Ta = Ta(t) is the air temperature.
The vertical boundaries x = 0, X are free from thehorizontal components of the velocity Vx and heat flow
Fx
Vx = 0, Fx = 0;

Fig. 3 Location of zero isotherm for various Tw(7)

In the considered geographical region, the meanannual temperature is positive, therefore, the zero
isotherm moves away from the upper surface andreaches the steady-state position in 10–15 years. As
seen, the variation of the temperature Tw inside theproduction-justified range has little influence on thethawing depth.

at the boundary z = Z
T = Ti .

(8)

The system of the differential equations (1)–(3)with the boundary conditions (4)–(8) and the initial data
T(0, x, z) = Ti , P(0, x, z) = 0 was solved using the finite
element method Zienkiewicz et al. (2013) and theoriginal codes (Nazarova, 1985; Nazarova et al., 1998).
The linear size of an element was 2 m, the time step was
6 h. The calculation involved the following model
parameters (Fig. 2): X = 600 m, Z = 100 m, z0 = 2 m,
ρ = 1020 kg/m3. The values of the coefficients ofconsolidationandtemperature
diffusioncorresponding to weakly permeable clays withmedium moisture content Goncharov (2002) arecompiled in the Table 1.

Fig. 4 Level curves of pressure (bar) 10 years after
filling of the tailing pond: a) undisturbed impervious
screen; b) impervious screen with rupture.

3

10th Asian Regional Conference of IAEG (2015)

Figure 4a illustrates distribution of the pressure P
at the time t=10 years when the length r (refer toFig. 2) of the horizontal boundary section A0A1 of theimpervious screen is 5 m : evidently, infiltrate reaches
the surface of the right-hand slope of the protective
embankment by now.
With this end in view, we analyze the influence ofthe value of r on the amount of infiltrate having
reached the outer perimeter of the embankment
L=A7A8A9A10A11 (refer to Fig. 2):t

Q (t )  Yvarious depth in these holes. The temperature
measurements are mainly required to monitor
propagation depth of the thawing zone whereas they
do not inform on the rupture of the impervious screen
since the temperature is almost the same at the bottomof the tailing pond.
Let us consider how the pressure field “responds”
the rupture of the impervious screen. At the time t = that the point x = xh, due to low-quality sealing
Ishchenko (2007), a hole 2c in length forms in theimpervious screen. In this case, in the section[xh–c, xh+c] of the upper boundary of thecomputational domain, the hydrostatic pressure ρgzh is
set so to be in accord with the water body depth zh at
the given point (Fig. 2). Evidently, the hole is anadditional source of fluid, and the infiltration toward
the daylight surface of the embankment intensifies(Fig. 4b) with closer location of the hole to theembankment surface (dashed colored lines in Fig. 5)
despite the fact that the water head is twice as little at
xh = 230 m than at xh = 170 m.
The colored lines in Fig. 6 show time variation of
the pressure Ph(t) =P(t, xm , zm) at c = 1 m; thecoordinates xm and zm of the observation points as wellas parameters th and xh of hole are specified in theinsets. For the comparison, the black lines show thepressure at the same points for the undamaged
impervious screen. It turns out that the pressure
sensors perceive the damage of the impervious screen
within the first months (depending on the distancefrom the damage position) to the accuracy of 10 kPa.
The slope α of the curve Ph(t) grows with closer
distance to the hole in the screen. The behavior α is theindication of the integrity of the impervious screen.

 V  n dL .
0 L

The varicolored lines in Fig. 5 show the change in
Q for different r as against Q = Q0 at the timet = 10 years when r = 0. It appears that a comparatively
slight growth of r essentially reduces the infiltrate
amount and, consequently, mitigates ecologicalhazard.

Fig. 5 Relative
embankment

volumeofinfiltrate

through

The numerical experiments have shown that Tw
has nearly no influence on the value of Q, while thedecrease in the mean annual temperature greatly
reduces the infiltrate amount due to the longer timewithin which the frozen layer occurs at the daylightsurface of the embankment.
3.

Search the location and timing ofimpervious screen rupture basedon piezometric measurement data

To monitor the state of the protective embankment,
vertical observation holes are drilled (see Figs. 1
and 2), and temperature and pressure are measured at

Fig. 6 Pressure at observation points for different
location xh of rupture in the impervious screen

410th Asian Regional Conference of IAEG (2015)

Now, we formulate the inverse problem: search
the time th and rupture point xh and in the impervious
screen by the data of piezometric measurements. Let
introduce an objective function: (th , xh )  I (th , xh ) / t* pa ,

(9)

where
t m  t*I (th , xh )  [ P(t , x , z , t , x )  P (t )] dt , (10)
2

m

m

h

h

*

tm

Fig. 7 Isolines of the objective function Ψ at theobservation point with the coordinates xm = 150 m,
zm = 32 m.

The function Ψ is a dimensionless functional of
discrepancy between the pressure P*(t) measured at
the observation point (xm , zm), that is the input data,
and the “theoretical” pressure P(t, xm, zm, th, xh) at thesame point, additional arguments of the function P
imply that this is the direct problem (1)–(8) solution atcertain values of th and xh. In (9), (10) tm is the start
time of the pressure growth at the point (xm , zm); t* isthe time interval set empirically. The numericalexperiments showed that t* grows as kS reduces; at thechosen model parameters, t* = 4–6 months.

It is found that Ψ is not always unimodal (Fig. 8);
for this reason, the objective functions are to beconsidered for various observation points.
The equivalence domain Um (yellow-colored in
Figs. 7 and 8) is a variation range of the objectivefunction arguments such that the relative change|ΔΨ|/Ψ is not higher than a certain small value (in ourcase, 0.005). The dimensions of Um , which condition
the accuracy ε of the inverse problem solution, dependon the noise level δ in the input data. The given
analysis assumes ε = 0.1, then the related value of δ is

The minimum point (th* , xh* ) of the objectivefunction furnishes the inverse problem solution. In
order to examine its resolvability, we synthesize inputdata:

P* (t )  [1   (t )]P0 (t ) ,where P0(t) - the exact solution of the direct
problem - is the pressure at the point (xm , zm) if theimpervious screen gets ruptured at the point x  xh0 at
the time t  th0 ; ξ is a random variable uniformly
approximately equal to 0.4.

distributed over the interval [-1,1]; δ is the relative
amplitude of superimposed multiplicative noise.
Figures 7 and 8 show the level lines of the objective

Fig. 8 Isolines of the objective function Ψ at theobservation point with the coordinates xm = 210 m,
zm = 30 m.

function Ψ for δ = 0.3, t* = 0.5 year, xh0 = 160 m,

The inverse problem solution ( xh* , t h* )  U , whereth0  4 years for different xm and zm .

U is the intersection of the equivalence domains Um of
the objective functions constructed for a number of5

10th Asian Regional Conference of IAEG (2015)

points where the pressure is measured. Then, with themoderate noise in the input data, the inverse problem
is unambiguously resolvable.
Based on the results of the numerical experiments,
it is expedient to use gradient methods Avriel (2003),
Nazarov et al. (2013) to search the minimum of Ψ if
the measurement points are sufficiently distant fromthe frozen zone since coefficients of thermal
diffusivity and filtration are discontinuous functions in
the vicinity of the zero isotherm and the objectivefunction Ψ becomes nondifferentiable as aconsequence.
4.

References

Avriel, M. (2003): Nonlinear Programming: Analysis
and Methods. Dover Publishing, 532 p.
Bejan, A. (2013): Convection Heat Transfer. Wiley,
4th Edition, 696 p.
Bense, V. and Person, M. (2008): Transienthydrodynamics in inter-cratonic basins during glacial
cycles, J. Geophys. Res., 113, F04005.
Buiskikh, A.A. and Zamoshch, M.N. (2010):
Prediction of thermal regime within a tailing dump
under permafrost, Journal of Mining Science,Vol. 46, No. 1, pp. 28-33.
Chzhan, R.V. (2002): Temperature condition andstability of low-head hydraulic power systems andground channels in permafrost zone. Yakutsk:
Institute of Permafrostology SB RAS, 207 p. (inRussian)Goncharov, S.A. (2002): Thermodynamics. Moscow:
Publishing House of Moscow State Mining
University, 440 p. (in Russian)Goy, L., Fabre, D. and Menard, G. (1996): Modelingof rock temperatures for deep alpine tunnel projects,
Rock Mechanics and Rock Engineering, Vol. 29,
No. 1, pp. 1-18.
Ishchenko, A.V. (2007): Infiltration protection andimpervious
efficiencyofwaterworks.
Rostov-on-Don: North Caucasian Research Center ofHigher School, 256 p. (in Russian)
Lolaev, A.B., Spesivtsev, A.V., Spesivtsev, V.V. et al.
(1997): Site investigation of tailing dam inpermafrost region, Geoenvironmental Engineering.
eds. R.N. Yong and H.R. Thomas. Thomas Telford,
London. 657 p.
Nazarov, L.A., Nazarova, L.A., Karchevskii, A.L.,
and Panov, A.V. (2013): Estimation of stresses anddeformation properties of rock masses based on thesolution of an inverse problem by the data of themeasurement of free surface displacements, Journalof Applied and Industrial Mathematics, Vol. 7, No. 2,
pp. 234-240.
Nazarova, L.A. (1985): Stress state of a
sloping-bedded rock mass around a working, Journalof Mining Science, Vol. 21, No. 2, pp. 132-135.
Nazarova, L.A., Nazarov, L.A., and Leontiev, A.V.
(1998): Three-dimensional geomechanical model of
the Tashtagol iron-ore deposit, Journal of Mining

Conclusions

The authors have designed the geomechanicalmodel describing the heat-and-mass transfer processes
in the vicinity of the protective embankment at thetailing pond of the Kumtor gold deposit, Kyrgyzstan,situated in the permafrost zone. The numericalexperiments have shown that: a) the zero isotherm
between the frozen and thawed rocks attains stablelocation 10–15 years after the embankment creationand the pond filling; b) the decrease in the meanannual temperature essentially reduces infiltrate
volume due to longer time when impermeable frozen
layer covers the embankment slope surface; c) thechange of the fluid temperature in the
production-justified range has almost no influence on
the depth of thawing and the amount of infiltrate in theembankment; d) a slight increase in the size of thehorizontal section of the impervious screen greatly
reduces the amount of infiltrate and, consequently, thepollution of the area.
The inverse problem on locating and timing arupture of the impervious screen is unambiguously
resolvable at an admissible accuracy for any real
thermal and filtration characteristics of soil using thedata of piezometric pressure measurements taken at anumber of points and given the moderate noise in the
input data.
Acknowledgements

The authors are grateful for the support provided
by the Russian Foundation for Basic Researches(Project No. 14-05-90116).

6

10th Asian Regional Conference of IAEG (2015)(2012): Increasing subsurface water storage indiscontinuous permafrost areas of the Lena Riverbasin, Eurasia, detected from GRACE, Geophys. Res.
Lett., 39, L09403.
Vozin, V.F. (Ed.) (1987): Environmental impact of ahydraulic power plant in the far north. Yakutsk. 88 p.
(in Russian)Yershov, E.D. and Williams, P.J. (2004): General
geocryology. Cambridge University Press, 608 p.
Yi, S., Arain, M.A. and Woo, M-K. (2006):
Modifications of a land surface scheme for improved
simulation of ground freeze-thaw in northern
environments, Geophys. Res. Lett., 33, L13501.
Zhang, Y., Chen, W., and Riseborough, D.W. (2008):
Disequilibrium response of permafrost thaw toclimate warming in Canada over 1850–2100,
Geophys. Res. Lett., 35, L02502.
Zienkiewicz, O.C., Taylor, R.L. and Zhu, J.Z. (2013):
The Finite Element Method: Its Basis andFundamentals. 7th Edition, Butterworth-Heinemann,
756 p.

Science, Vol. 34, No. 3, pp. 209-216.
Nepomnyashchaya, T.A. (2007): Macro-zoobenthosof the Kumtor river under anthropogenic impact,
KRSU Bulletin, Vol. 7, No. 1, pp. 111–113. (inRussian)Osmonbetova, D.K. (2011): Modern ecology in the
Naryn river upstream zone and predictive estimatesin view of the gold deposit development, KRSU
Bulletin, Vol. 11, No. 3, pp. 50-53. (in Russian)Rawlins, M.A., Ye, H., Yang, D. et al. (2009):
Divergence in seasonal hydrology across northernEurasia: Emerging trends and water cycle linkages,
J. Geophys. Res., 114, D18119.
Riseborough, D., Shiklomanov, N., Etzelmuller, B. etal. (2008): Recent advances in permafrost modeling,
Permafrost and Periglacial Processes, 19,
pp. 137-156.
RF Construction norms and regulations 2.01.01-82(1996): Construction climatology and geophysics,
Ministry of Construction of Russia. 140 p. (inRussian)
Velicogna, I., Tong, J., Zhang, T. Kimball, J.S.

7