Numerical Modelling of Tailings Dam Thermal-Seepage Regime Considering Phase Transitions

Statement of the Problem. The article describes the problem of combined thermal-seepage regime for earth dams and those operated in the permafrost conditions. This problem can be solved using the finite elements method based on the local variational formulation. Results. A thermal-seepage regime numerical model has been developed for the “dam-foundation” system in terms of the tailings dam.The effect of heat-and-mass transfer and liquid phase transition in soil interstices on the dam state is estimated.The study with subsequent consideration of these factors has been undertaken. Conclusions. The results of studying the temperaturefiltration conditions of the structure based on the factors of heat-and-mass transfer and liquid phase transition have shown that the calculation results comply with the field data. Ignoring these factors or one of them distorts the real situation of the dam thermalseepage conditions.


Introduction
Today, intensive development of the northern regions of Siberia and the Far East of Russia is associated with the exploitation of their rich natural resources.This development requires operation of already built hydraulic engineering structures of the power and mining industries, water supply facilities, and construction of new ones.Herewith, we need to deal with the problems of construction in the permafrost conditions [1][2][3].Earth embankments and dams are the most frequent hydraulic engineering structures used in the permafrost conditions [1,3].Trouble-free operation of these structures depends on compliance with the required thermalseepage regime of the dam-foundation system [2][3][4].
Any disturbance of this mode may cause emergency and even structure collapse [4][5][6][7].According to experts, about 50% of 3rd and 4th category earth low-head dams built in the severe climatic conditions is damaged during the first three years of their operation due to thermal-seepage regime disturbance [4].
When designing the earth hydraulic engineering structures in the permafrost conditions, the extremely important task is the reliable determination of the thermal-seepage regime of the dam-foundation system.In most cases, the geotechnical setting in the construction area is nonhomogeneous geology and composed of frozen and thawed layers stratification.In order to determine the exact seepage pattern at the hydraulic engineering structure foundation, it is often necessary to solve the spatial three-dimensional problem [8][9][10][11].In some cases, the seepage heating effect may result in negative consequences, such as increasing soil permeability and seepage losses in a reservoir, reduction of bearing capacity of frozen rocks, and seepage failure [1].
Therefore, for the facilities located in the area of permafrost occurrence, handling the problem of the combined thermal-seepage regime is actual based on the phenomena of heat-and-mass transfer and phase transitions, whose results are to be considered in the design.
Both in Russia and abroad there are numerous scientific studies devoted to the problems for calculations of seepage, thermal, combined thermal-seepage regimes of earth dams and their foundations [12][13][14][15][16].In recent years, a mathematic modelling is widely used for these purposes [5,13,[16][17][18].However, most software systems including industrial ones handle the seepage and temperature problems with no interaction of two processes [8,[19][20][21][22][23].Additionally, the effect of liquid phase transitions in various soils on thermal-seepage regime formation is not entirely clear.Thus, the need for further scientific studies in this sphere including calculation method improvement is quite obvious.The purpose of this work is to improve the method for calculating the thermalseepage regime of the dam-foundation system based on heatand-mass transfer and phase transitions in soil interstices.

Method
For the classical pattern, the temperature problem (Stephen) with no heat-and-mass transfer in the seepage flow can be solved by the differential equation [1]: where  = (, , , ) is a temperature function;  is a time;   ,   , and   are temperature conductivity coefficients in -, -, and -directions.
The set of equations describing dam thermal regime considering the temperature of seepage water is known from the heat-and-mass transfer theory: where  = (, , , ) is a foundation or dam material temperature; (, , , ) is a seepage water temperature;  is a time;   ,   , and   are conductivity coefficients of dam body or foundation material temperature in -, -, and -directions; V  , V  , and V  are seepage rates in -, -, and -directions;  V and   V are the coefficients that characterize the heat exchange between the material of the dam or foundation and filtering water in its pores.These coefficients are determined from the following dependencies: where  V is a volume heat-exchange rate specifying heat exchange between the dam or foundation material and ambient water;   and   are volumetric specific heat of water and solid body;   and   are water and solid body density.The solution of the set of equations (2) describes change of seepage water temperature after passing through the dam body or foundation material.Heat exchange in each space point is featured by a volumetric heat-exchange rate.The calculation does not use the coefficient of water temperature conductivity due to insignificant conductive heat transfer.
If the set of equations ( 2) is considered with regard to soil porosity, then the set changes to [5] where  is a soil porosity.
If the problem contains a temperature of dam body material or foundation soil equal to a seepage water temperature, then the combined problem is the solution of Fourier-Kirchhoff equation [5]: where For soils with minimum porosity, (5) changes to Search for solving the combined thermal-seepage problem comes to solving differential equations (( 5), ( 7)) for the area with seepage liquid flowing and (1) with no seepage.
The differential equations solution comes to minimization of functional  [5].The standard form of (7) based on boundary conditions and phase transition is given below: where scalar product of seepage velocity vectors in -, -, and directions and temperature gradient; Ω 1 , Ω 2 are the surfaces with the designated 2nd and 3rd category boundary conditions;  phase is the temperature of phase transformations.
The proposed method for calculating the combined thermal-seepage equation based on heat-and-mass transfer and phase transition has been implemented by the authors in the FILTR_FAZ software system.This software allows solving both classic filtration tasks and heat conduction problems, as well as obtaining complex combined solutions.To date, there are also universal industrial computing systems (such as ANSYS [20], ANSYS CFD, MIDAS CFX, Feflow, and Plaxis Flow), but at the same time, their use is only possible to search for the simplest classical solutions.Without taking into account the complex physical processes taking place in permafrost, using FILTR_FAZ the authors have studied the thermal-seepage regime of hydraulic engineering structures in the permafrost (e.g., [5]).Some results of these studies are given below.

Problem Statement and Content.
The facility is a thawed type hydraulic fill dam in terms of its design.The dam upstream fill is erected by aggregation, while the downstream one is filled with soil from open pits.The downstream face is fastened by means of mined rock with an effective thickness of 3.5 m.The medium fine gravel entrenchment is used as an anti-seepage cutoff.
The structure is erected in the severe climatic conditions.The average winter temperature is minus 24.8 ∘ b and the summer temperature is plus 15.6 ∘ b.The recorded absolute minimum is minus 44 ∘ b, the absolute anomalous maximum is plus 33 ∘ b, and the snow depth is over 1 m.The foundation permafrost thickness is over 100 m.Yearly average air temperature (according to the field studies) is within the range of minus 1.0 to minus 1.5 ∘ b.
The embank height is 72.00 m with an edge elevation of 360.00 m.The dam is expected to be topped by 2032 up to an edge elevation of 370.00 m.Therefore, the structure height will reach 82 m.This facility belongs to the 1st importance class.
The 3D simulation mathematic model has been developed to determine the combined thermal-seepage regime of the facility (see Figure 1).This model describes an earth dam, retaining prism, drainage system, lower manoeuvring basin, and pond beach.
The combined thermal-seepage regime was calculated by an improved procedure with the use of the FILTRT_FAZ and FILTR software system [5].The phase transitions of interstitial water and heat-and-mass transfer in the direction of seepage flow were taken into account in the calculations.Multipurpose industrial software ANSYS APDL for mathematic modelling was used to analyze the seepage regime in process.
The purpose of the study was to determine the combined thermal-seepage regime in the 3D pattern based on the phenomena of heat-and-mass transfer and phase transitions within 10 years of the date of the facility erection up to a design elevation of 370.00 m.
To analyze the structure condition, the following problems were solved: seepage problem for a design area of 2.5 km 2 ; temperature problem in its classical pattern; combined thermal-seepage problem with due regard for heatand-mass transfer by seepage flow; combined thermalseepage problem based on heat-and-mass transfer by seepage flow and phase transitions of interstitial water.

Seepage Problem Solution.
The calibration tests are carried out at the initial stage to determine the position of depression curve in the model similar to piezometers reading.The calculations were used to determine P  seepage coefficient for the retaining prism and paving rock-fill blanket.The study results have shown that the retaining prism is most likely mudding at the upstream side, and, therefore, two P  seepage coefficients were used in the calculations; that is, P  for the mudding section is 50 m/day and P  for other sections of the rock-fill blanket is 100 m/day.The seepage calculation results in dam section 1-1 are shown in Figure 2.
The seepage regime calculations based on the classical pattern have shown that the primary seepage flow is directed through the left-bank abutment (see Figure 3), which is confirmed by the field data.
In order to more accurately determine the depression surface position, it was identified iteratively with control of the finite-element grid within its expected location.
The seepage regime was corrected at every stage by time based on changes in the thermal regime of the previous stage.The lower values of P  were obtained for the elements, whose temperature values of the thermal equation were within the range of minus 1.0 ∘ b and below and set by the absolute confining layer with temperatures of minus 0.5 to minus 1.0 ∘ b.Then the correction calculation was performed.The resulting seepage velocities at every stage were used to determine the structure thermal-seepage regime.

The Thermal Problem Solution Based on Heat-and-Mass
Transfer and Phase Transitions.The initial and boundary conditions were set to solve the thermal problem.The initial condition (temperature distribution in the soil body at the calculation) was set according to the field data as of November of the current year.The condition of convective heat exchange (the 3rd category boundary condition) was set according to soil-air contact.The mean monthly air temperatures are given in Table 1.The 1st category boundary condition,   =  pr (where  pr is minus 1.5 ∘ b and  pr is the temperature of permafrost), was set in the permafrost location.The absolute thermal insulation was set according to the boundaries of design area (the special pattern of the 2nd category boundary condition).The 1st category boundary condition,   =  water , where  water is a water temperature in every design month, was specified according to the soil-air contact.
The thermal-seepage regime was successively calculated for every month.
The thermal problem was solved in several stages for 10 years' period of structure operation.The solution of absolute thermal problem was obtained at the first stage without considering heat-and-mass transfer and phase transitions.
The temperature patterns for dam I-I section in November for 1, 3, and 10 years of operation are shown in Figure 4.
In this case, a frozen core with temperatures of 0 to minus 0.5 ∘ b (higher than in the dam body) is formed in the centre of the thawed dam due to the boundary conditions being set according to the field data.The downstream face freezes to the permafrost rock top forming a negative temperatures field at the downstream side.
According to the temperature distribution in the soil body, it may be concluded that the results obtained from the classical pattern of the thermal problem do not comply with the field data.Thermometers located on the crest record positive soil temperatures ranging from 0 to plus 0.5 ∘ b.The seepage problem solution for this option showed the seepage flow occurring on the surface of the downstream face only, that is, in the area of seasonal freezing/defrosting.Also this fact does not comply with the true pattern of the seepage regime.
The solution of the combined thermal-seepage problem was considered at the second stage based on the heatand-mass transfer phenomena.The temperature distribution pattern partly corresponding to the field data was obtained according to the calculation results (see Figure 5).The positive temperatures field ranging from 0 to plus 0.3 ∘ b is formed in the dam centre, but the temperature field between the thawed core and retaining prism is nearly unchanged.This is provided by severe exposure to winter negative temperatures that causes freezing to the soil body at the downstream face.Therefore, the seepage flow remains trapped in the dam centre, where "a pocket" is formed.In case of soil thawing in summer along the downstream face, high rate seepage takes place and it has no impact on the negative temperatures field smoothly flowing it around.
The defrosting effect due to heat-and-mass transfer occurs in the seasonal thawing area only.In comparison to the previous option, this area increases in size, but it completely froze in winter.
At the final stage of calculation, the combined thermalseepage problem has been solved based on phase transitions of cohesive water in the soil (see Figure 6).
In the first years, while considering phase transitions, the soil body will be defrosted in the dam centre (as confirmed by instruments reading).A thawed area with passing seepage flow is formed beneath the layer of seasonally frozen soil.The thickness of the seasonally frozen area is 3.5-3.7 m that complies with climatological data of the region.
The frozen area from the downstream face formed at the initial stage gradually thaws out under the influence of heatand-mass transfer and phase transitions of cohesive water.Since the foundation soils are characterized by the negative temperatures ranging from minus 1.0 to minus 1.5 ∘ b, the heat carried by the seepage flow is sufficient for soil thawing.Based on the fact that the structure was designed by the thawed type, this operation mode corresponds to the design one.
Retaining prism permafrost thawing is observed at the downstream face by the 10th year of operation causing additional setting.The most extensive temperature change is observed in the first 6-7 years of dam operation, and then the temperature field tends to be stabilized.However, some local areas of the negative temperatures in the thawed body remain by the 10th year of dam operation.
The key feature defined during the design studies is fast degradation of the permafrost in the way close to the soil phase transitions temperature; it can be explained by the concept of incomplete water freezing in dispersive frozen soils [24,25].At any temperature rise (including the negative values field), frozen soils thaw out.This phenomenon caused accidents of many structures built in the permafrost conditions.
The results of problem solution as compared to the field data obtained in a dam temperature well for 2015 are given in Table 2.
Difference between obtained data and instrument readings is 0.5 ∘ b maximum that is less than 0.1% in terms of the annual variations amplitude.When considering phase transitions, the soil condition in the dam body may be exactly determined in contrast to other solutions.
The following conclusions are made based on the results of the calculation studies: (1) The improved calculation procedure and software system FILTRT_FAZ created by the authors allow handling the combined thermal-seepage problems based on the phenomena of heat-and-mass transfer and phase transitions.
(2) The results of studying the temperature-filtration state of the structure based on the factors of heat-and-mass transfer and liquid phase transition have shown that the calculation results comply with the field data.Ignoring these factors or one of them distorts the real situation of the dam temperature-filtration conditions.

Figure 1 :
Figure 1: Solid approximant of the earth tailings dam.

Figure 2 :Figure 3 :
Figure 2: Calibration test results for the dam seepage regime at two P  values of the rock-fill blanket: 50 m/day and 100 m/day.

Figure 4 :
Figure 4: Temperature field in dam I-I section as of November according to the classical pattern of the thermal problem solution: (1) in 1 year of operation; (b) in 3 years; (c) in 10 years.

Table 1 :
Mean monthly temperatures of ambient air and water.