Numerical analysis of tailing dam with calibration based on genetic algorithm and geotechnical monitoring data

: The article presents a method of calibration of material parameters of a numerical model based on a genetic algorithm, which allows to match the calculation results with measurements from the geotechnical monitoring network. This method can be used for the maintenance of objects managed by the observation method, which requires continuous monitoring and design alterations. The correctness of the calibration method has been verified on the basis of artificially generated data in order to eliminate inaccuracies related to approximations resulting from the numerical model generation. Using the example of the tailing dam model the quality of prediction of the selected measurement points was verified. Moreover, changes of factor of safety values, which is an important indicator for designing this type of construction, were analyzed. It was decided to exploit the case of dam of reservoir, which is under continuous construction, that is dam height is increasing constantly, because in this situation the use of the observation method is relevant.


Introduction
Geotechnical monitoring is closely linked to the observational method, which was described by Peck (1969). It assumes constant design changes resulting from information obtained from measurements. Nowadays, data collection and processing are becoming easier even with large datasets. Design modifications require calculations that are less and less time-consuming due to the development of computers and implementation of numerical methods. Moreover, numerical methods allow to calculate complex structures as a whole without the need to separate smaller calculation schemes. These methods also allow to take into account poromechanical coupling, which is extremely important for challenging geotechnical issues, especially in case of dams.
Despite the constant development of computer programs, some aspect of FEM modelling remain in the designer's possession. The selection of material parameters is a useful example. It was stated by Gioda and Sakuri (1987) that back analysis is a practical tool for interpreting geotechnical monitoring data and for determining the material parameters. The same concept is used by researchers (e.g. Gioda and Sakuri, 1987;Zentar et al., 2001): the calibration incorporates parameter changes to fit the results from the model to the measurement data.
The application of the genetic algorithm for this purpose has been successfully presented in many publications. The results obtained by Vahdati, P. et al. (2013Vahdati, P. et al. ( , 2014 show that this type of optimization works well for the selected parameters on the basis of data gathered from one gauge placed on a real structure for elastic-perfectly plastic Mohr-Coulomb model and for more sophisticated Hardening-Soil model. Other publications have proved that genetic algorithm can be applied to optimize material parameters in laboratory test (Rokonuzzaman and Sakai, 2010) and to predict the behavior of soil in future, (Papon et al., 2011) based on fitting numerical results to data from pressuremeter tests and resonant column testing. However, most construction maintained with observational method are large and data is collected from several sensors. Moreover, in such a situation, boundary condition can be time-dependent and data used to fitting should be treated as a function of time, which was presented by Gioda and Sakuri (1987).
According to Whitley (1998) a genetic algorithm is any optimization model that uses selection and recombination to create new points in the searched space. It should be noted, however, that most models are based 2 Szczepan Grosel on the so-called classic genetic algorithm (Holland, 1975). In this paper, it was decided to use an algorithm based on selection by ranking. The publication (Whitley and Starkweather, 1990) shows that such an approach may be more effective than the classical one.
In this work the methodology of applied optimization algorithm and the result of calibration carried out for five different materials, which yields a total of 25 parameters, based on the data from 10 sensors is presented. Synthetic data was used to estimate the correctness of the algorithm. Also prediction, that is calculation concerning the future behavior of the model, are considered to check the suitability of using the optimization algorithm as a complement to the observational method. It was decided to use the numerical model of tailing dam as an example. Such geotechnical constructions are often complex with the crucial impact of filtration. Some of them are also maintained with observational method, for example OUOW Żelazny Most in Poland (Jamiolkowski, 2014).

Methodology of optimization algorithm
The calibration method is based on a genetic algorithm, that is it uses recombination to create new solutions. When optimizing the parameters of the numerical model, it is necessary to perform calculations in order to select the best solutions. These calculations take much longer than the rest of the algorithm. In addition to the aforementioned advantages of the ranking approach, it also minimize the calculation time, because less numerical calculations need to be performed. The description of the algorithm is divided into subsequent steps, explained in the following subsections. The algorithm has been implemented using Python environment (Van Rossum and Drake, 2009; Oliphant, 2006;McKinney, 2010;Hunter, 2007). Due to the fact that this type of calibration will not always find the optimum, it was decided that its parameters (mentioned below) are not fixed, but should be selected to get the best results for the specific task. Every calibration user has to set these parameters according to his experience.

Measurements data
Prior to performing the optimization, it is required to specify the data to which numerical model should be fitted. In this work, it was assumed that measurement data are functions of time, for example changes of displacement in time (data from benchmarks) or changes of piezometric head in time (data from piezometers).

Numerical (FEM) model
The numerical model should be created using all possible construction information, such as geological layers arrangement, history of loading or geotechnical tests. The size of the model largely determines the time needed for optimization. Selection of nodes or elements that correspond to actual sensors is necessary to perform the algorithm. The results needed from the model are only appropriate values (displacement, pressure) of the selected nodes or elements in text form. It is obligatory to use FEM software that allow to create such an output automatically, that is without running the postprocessor.

Individuals encoding
The genetic algorithm is based on population that consist of individuals. The individual represents the parameters of all optimized materials. As the Mohr-Coulomb model was adopted to numerical calculation, parameters for each material include: angle of internal friction, cohesion, Young's modulus, Poisson's ratio and hydraulic conductivity. Different sets of parameters can be selected for each material. The individual is a sequence of numbers, which determines all the parameters according to the following scheme: where: -angle of internal friction, -cohesion, -Young modulus, ν -Poisso conductivity, -number of optimized materials, and the subscript refers to mate hole dataset according to following formula: where: -number of gauges in eme: where: ∨, ∧ -crossover point, 1, 0, , -individuals features. ern: where: 1, 2, 3 -three best ranked individuals. ern: , , , , , , , … , , scheme: where: -angle of internal friction, -cohesion, -Young modulus, ν -Poisson ratio, -hydraulic conductivity, -number of optimized materials, and the subscript refers to material numbe hole dataset according to following formula: where: -number of gauges in eme: where: ∨, ∧ -crossover point, 1, 0, , -individuals features. ern: where: 1, 2, 3 -three best ranked individuals.

ern:
where: φ -angle of internal friction, c -cohesion, E -Young modulus, ν -Poisson ratio, k -hydraulic conductivity, n -number of optimized materials, and the subscript refers to material number.

Initial population
The initial population of individuals is generated randomly. However, parameters are drawn from ranges specified by the designer. First reason for such an approach is the fact that some parameters cannot exceed certain values, for example, the Poisson's ratio. Secondly, these ranges can be determined on the basis of geotechnical documentation, which can lead to a result consistent with engineering experience. Additionally, the calibration time can be reduced.

Evaluation
To perform the evaluation, numerical calculation must be carried out. The same task is calculated, but the optimized materials parameters are assigned according to the individuals. Evaluation of an individual involves calculation error defined as the difference between measurement data and numerical output. Since measurement data can have different physical sense and therefore different units and values, data should be normalized or relative error can be used. The latter was implemented, the error is based on mean relative error estimated on the whole dataset according to the following formula: gle of internal friction, -cohesion, -Young modulus, ν -Poisson ratio, -hydraulic -number of optimized materials, and the subscript refers to material numbe ccording to following formula: , mber of gauges in eme: crossover point, 1, 0, , -individuals features.
-three best ranked individuals.

cement amo
where: r -number of gauges in calibration, n i -number of measurements in i-th gauge, x n -measurement value from numerical analysis, x m -measurement value from monitoring.

Selecting individuals for reproduction
Three individuals with the lowest error values are selected for reproduction. In the method (Whitley and Starkweather, 1990) only two best individuals had been selected, however, it was decided to increase this number due to the computing capabilities of computers, which allow to calculate several FEM models at the same time, while maintaining the crossover of the best individuals.

Crossover
Crossover is the exchange of information between two individuals. The crossover point is selected, in relation to which the sequence of numbers is divided. Then the parts of individuals are swapped as shown in the following scheme: -angle of internal friction, -cohesion, -Young modulus, ν -Poisson ratio, -hydraulic ity, -number of optimized materials, and the subscript refers to material numbe et according to following formula: number of gauges in eme: ∧ -crossover point, 1, 0, , -individuals features.
As a result, two new individuals are created, but only one is randomly selected for the next population. Pairing of individuals for crossover occurs according to the following pattern: where: -number of gauges in eme: where: ∨, ∧ -crossover point, 1, 0, , -individuals features. ern: where: 1, 2, 3 -three best ranked individuals. ern: imposes replacement amo where: 1, 2, 3 -three best ranked individuals.
In each generation, i.e. that is every next population, the crossover point is randomly selected for each crossover. It should be emphasized, that the way parameters are written as individuals is important, although the parameters are independent, several of them describe one material. Already presented encoding scheme enforces the exchange of parameters in all materials, however, the following pattern: where: -angle of internal friction, -cohesion, -Young modulus, ν -Poisson ratio, conductivity, -number of optimized materials, and the subscript refers to material nu hole dataset according to following formula: where: -number of gauges in eme: where: ∨, ∧ -crossover point, 1, 0, , -individuals features. ern: where: 1, 2, 3 -three best ranked individuals. ern: imposes replacement amo imposes replacement among the materials. Since it is highly time-consuming to establish which method deliver the best results selecting the encoding system is treated as an algorithm parameter.

Mutation
Mutation introduces small changes to the newly created individuals. Maximum values of changes are defined for each type of parameter, for example for all cohesion values mutation parameter is the same. Mutation can be conducted as addition or multiplication. The second option is useful for hydraulic conductivity, as its values can differ by several orders of magnitude. Parameters which resulted from mutation cannot exceed the ranges defined for initial population. A certain number of parameters (defined as algorithm parameter) are randomly selected for mutation. Each parameter can be increased (addition or multiplication by a mutation parameter), decreased (subtraction or division) or remain unchanged due to the random choice.

Next population
Each subsequent generation consists of individuals from the previous generation without the three worst ones; in their place three individuals created by crossbreeding and mutation are added. In the first generation, FEM analysis must be carried out for all the individuals, while in each subsequent generation, for only three new ones.

Synthetic measurement data, numerical calculation and algorithm's parameters.
In this chapter the results of calibration carried out for synthetic data are presented. A numerical model with some material parameters was created and calculated. Selected nodes were treated as sensors and displacement and piezometric head values over the time were recorded. The data prepared in this way were noised by adding values from normal distribution with mean equal to 0 and standard deviation adapted to values of the data. Figure 1 presents an example of the noised synthetic data. The same model was used both to create the data and to perform calibration to ensure that a solution of optimization exists. Moreover, the impact of model inaccuracies was eliminated. The calculation problem concerned embankment of flotation reservoir. The height of embankment increased in time as the water and sedimentation particles filled the reservoir. Numerical model was created using FEM software ZSoil (Zimmermann et al. 2016) with assumption of plain strain. The first order EAS (enhanced assumed strain) (Simo and Rifai, 1990) quadrilateral elements were used. After 205 computational step the height of the dam was equal to 65 m, while the depth of subsoil was 265 m. In the figures 2 to 6 subsequent stages of construction history, geometry, mesh and model's dimensions are presented.
On both right and left edges the horizontal displacement was blocked, while on the bottom edge horizontal and vertical displacement was fixed. On the left edge, water pressure was raised to simulate reservoir filling, and the right edge water pressure had constant value throughout the whole calculation. As it was already mentioned elastic-perfectly plastic Mohr-Coulomb with unassociated plastic flow rule was assumed. Coupled analysis, i.e. that is consolidation, with adoption of van Genuchten's model of unsaturated zone (Van Genuchten, 1980) was selected to calculate the distribution of pore pressure and deformation of the structure. The model consisted of almost 12,000 elements and the calculation took less than 30 minutes (on 3.2 GHz CPU and 32GB RAM computer). Nodes selected as sensors are shown in Figure  7. Sensors with letter U are benchmarks, which gather information about vertical and horizontal displacement,       Figures 2 to 6. It was decided that 10 parameters could be mutated. Ranges used for the assessment of initial population as well as mutation parameters are presented in Table 1. Such a calibration took about 10 hours, of which only about 5 minutes were needed for the operations related to the performing genetic algorithm (not FEM analysis).

Results
The calibration was performed 5 times for different random initial populations, yet for the same initial ranges and mutation values. Parameter values selected to create synthetic data and optimized material parameters are presented in Table 1. Optimizations are described with Roman numbers I-V. It has to be mentioned that exact parameters have not been reproduced in any of the 5 optimizations. However good approximation was achieved for materials 2B, 3C and 5E. Surely, it is connected to role the played by each of these material in the model. Layer 2B corresponds to about 80% of elements in subsoil, layer 3C is the weak layer that enforces the shape of potential failure zone and layer 5E is the stiffest one. The hydraulic conductivity values have been optimized satisfactorily. Probably narrowing the ranges of Poisson ratio could   improve results, as these ranges were the same for all materials and relatively wide. Figure 8 shows the change of error depending on the generation number, that is how the error in each of the five calibrations was reduced. It can be observed that in all five cases error decreased, which means that the process gradually improves the initial individuals, which were created randomly. However, there is a possibility that there will be no significant improvement during the optimization process if the first population contains very good individual, similar to calibration IV. From a practical point of view, such a situation cannot be treated as proof of algorithm failure because the parameters provide a reasonable fitting to monitoring data. Figures 9 to 13 show the fitting diagrams for all the sensors. Data from measurement is marked with black color and in these plots they are shown without noise. On the X-axis the time as number of computational step is presented and on the Y-axis sensor quantity: displacement or piezometric head. It is difficult to decide which calibration (from I to V) is the best, but all of them should be acceptable in real case. It can be observed that for different sensors different sets of parameters correspond to the best fitting obtained.

Predictions
The fitting diagrams presented in the previous chapter show that the results obtained from all 5 parameter sets are similarly matched to the measurement data. Although the fits are not perfect, at the engineering level they are at least satisfactory. In other words, if the actual parameters on an existing object were not known, all 5 sets would be accepted as sufficient matches and all discrepancies would be treated as inaccuracies resulting from the fact that the numerical model is only an approximation of reality. However, from the designer's point of view, displacement values or piezometric head alone are not sufficient to assess the safety of a construction. The factor of safety is most often used for this purpose.
Therefore, factors of safety were calculated using strength-reduction method (Cała and Flisiak, 2003) with reduction of both strength parameters, i.e. cohesion and angle of friction. Then FOS values were compared during the history of the dam's erection, at the end of the period in which the measurement were available, and also in the future, assuming a continuous dam increase of 0.4 m per 1 calculation step, which gave a total increment of dam height equal to 18 m. Further increase would cause the failure surface to exceed the model area or to be too close to its boundaries. Figures 14 and 15 show the geometry of the model at the final computation step and Figure 16 depicts the comparison of FOS values.
Despite different FOS values for different sets of parameters, the failure surfaces are similar. The deeper layer of 3C determines failure surface, regardless of the set of parameters. In the Figures 17-20, evolution of the failure surface due to increasing dam height is shown for dataset IV; however, the shape of slip surface is valid for all optimized datasets and case with actual parameters.
In addition, predictions of the values that were monitored are also shown. The Graphs 21-25 show these  predictions -on the horizontal axis there is time (in the form of calculation steps) and on the vertical axis is the measured value (displacement or piezometric height). A continuous line indicates the fit to the existing measurements and a dashed prediction.
Both the graph of the factor of safety and the predictions of the measured values show that, despite a good fit, the predicted values may differ significantly from the actual values. What is worth emphasizing, however, is that the set of parameters which had the smallest error during optimization shows predictions, on the basis of which (while maintaining appropriate safety factors) design modifications assumed by the observation method can be prepared.

Summary
This article shows that the calibration of material parameters of a numerical model for monitoring data using the genetic algorithm-based optimization method is also possible for more complex geotechnical issues and that it also works with a larger number of sensors that collects different data. An outline of the implementation of the optimization method and its performance on the example of artificially generated data is presented. Although the results of the matching were satisfactory, the optimization algorithm may need to be changed for specific cases. It is also obvious that the results of the matching depend on the correct creation of a numerical model.
The second part presents the differences between predictions created by calculating numerical models using optimized parameters and simulating the numerical           Numerical analysis of tailing dam with calibration based on genetic algorithm and geotechnical ... 13 model used for data preparation. The most important conclusion to be drawn is that not every good match guarantees a prediction corresponding to reality. The best prediction was obtained for a set of parameters, which gives the smallest error during calibration, so it is recommended to carry out several calibration processes and select the one with the best fit. Perhaps a larger number of generation could provide better results, however it was decided to shorten the time of optimization. It should be stressed that the task of calibration of parameters is the inverse problem. The most common way to improve the outcome of the inverse problem is regularization. Another method in the case presented in the article is to better define the ranges from which the initial population is drawn and modify the mutation so that correlations between individual parameters are taken into account. Additionally, in the case of analyses of future behavior of geotechnical objects, a statistical approach can be applied and conclusions can be formulated based on a group of potential solutions. When discussing the problem from the perspective of inverse problem, the appropriate choice of sensor position should be considered. It would be worth examining whether there exists such a group of sensors that would allow to obtain a unique solution. It will certainly be the subject of the author's future research, especially as it would facilitate the design of monitoring networks and maintenance of geotechnical objects using the observation method.
However, it should also be noted that the predictions for the measurements of the sensors taken into account are the same for all the cases for a sufficiently short prediction period. This is consistent with the observation method, which requires not only continuous monitoring but also verifying the assumptions made and modifying the project. A much more serious danger is evident when comparing the factor of safety values. Large discrepancies already exist during data matching (about step 200, i.e. before prediction). As with prediction, the smaller the  error during calibration, the closer the values of FOS were to reality. Only the set of parameters for which the minimum error was achieved reproduces FOS in a designsuitable manner. Although all the predicted FOS values are smaller than the actual, that is safe; this cannot be regarded as a rule because of a small test group. The failure surfaces were correctly reconstructed.