A comprehensive approach to the optimization of design solutions for dry anti-flood reservoir dams

A 3D terrain and subsoil model is created on the basis of design documentation, i.e. hypsometric map (of terrain) and results of geotechnical investigation, e.g. boreholes or soundings. Based on the contour map, a three-dimensional terrain mesh model is generated using appropriate CAD tools. It is then imported into the FE software as the top boundary of the 3D computational domain. The domain should cover at least the area of the dam extended by its total width in every direction. The depth of the subsoil in the model should be at least several times greater than the height of the designed dam. Along with the division of the domain into finite elements, the geotechnical layer system is reconstructed. The reconstruction consists in assigning the appropriate material and its parameters to each element. It is proposed to be carried out using the kriging technique. In general, it consists in minimizing the error variance of the interpolated field assuming the so-called semi-variogram function, in particular the correlation length (or radius) (Oliver & Webster, 1990). Within this approach, the elevation of layer bottom/ceiling is treated as a random variable (Krige, 1951; Oliver & Webster, 1990; Webster & Oliver, 1993). Its values are given in the boreholes and need to be interpolated in the entire domain. The use of the kriging technique makes it possible to create a probable three-dimensional arrangement of soil layers based on the data from limited number of boreholes or geological profiles.

The optimization of the extent of the anti-filtration barrier is carried out by a proper parametric analysis. This includes the depth of the barrier and the length of its anchorage in the abutments. The minimum values of these geometrical features need to be determined so that to provide the fulfilment of a few conditions described below. The main consideration in this context is the limit state HYD related to unfavourable phenomena resulting from water flow in the soil such as seepage erosion or soil boiling. To avoid possible occurrence of internal erosion (suffosion), the so-called critical filtration gradients must not be exceeded. Furthermore, the possibility of soil boiling, being a result of filtration stability loss, is verified during the numerical analysis of the overall dam stability conducted using the shear strength reduction (SSR) method (see step (3) described in subsection 2.3). This, however, requires the hydromechanical coupling to be taken into account at least in terms of the effective stress principle. The same applies to the UPL limit state. Moreover, regarding the assessment of filtration phenomena, the amount of water flow under and around the anti-filtration barrier in the abutments area is determined. The total discharge of water through the subsoil must be limited to the value specified in the design assumptions.

In verifying the HYD limit state condition due to the possibility of internal erosion caused by suffosion, it should be proved that: (1) γi⋅i≤icr {\gamma _i} \cdot i \le {i_{cr}} where i is the hydraulic gradient, i_cr is the critical values of hydraulic gradient for a given soil, and γ_i is the safety factor, which depends on the class of structure (determined according to national regulations and separate standards).

In order to verify condition presented above, initial-boundary value problem concerning transient flow in partially saturated medium should be solved. In practice, the van Genuchten model of soil retention curve (Van Genuchten, 1980) and the Richards hypothesis of gas phase continuity in the aeration zone are most often adopted (Truty et al., 2020). The distribution of the hydraulic gradient i is determined according to Darcy’s law (Scheidegger, 2020) based on the fluid velocity distributions obtained from numerical calculations. Within the saturation zone, the following formula applies: (2) i=vk, i = {v \over k}, where k is the hydraulic conductivity of soil and v is the fluid velocity.

Since the critical values of the hydraulic gradient may be different for each layer, it is difficult to evaluate the gradients directly from condition expressed by Eq. (1). For this reason, it is proposed to transform the condition (1) to the following form: (3) Fi=γi⋅iicr≤1, {F_i} = {{{\gamma _i} \cdot i} \over {{i_{cr}}}} \le 1, where the value of the parameter F_i should be interpreted as the multiplicity of the critical gradient. Then, F_i values can be interpreted in a global way, i.e. in the whole model: F_i ≤ 1.0 means that the condition of permissible hydraulic gradient values is met, while F_i > 1.0 indicates the substrate areas where internal erosion could potentially occur. The critical values of hydraulic gradients in different soil types are determined most often from empirical relationships broadly described in the literature. Among the most widespread, there are the formulas of Sichardt (1928) and Abramov (1952), with the latter being more conservative.

Hydraulic gradients and water flow intensity under and around the barrier can be estimated based on the solution of uncoupled problem. On the other hand, it is required by Eurocode 7 to take the variation in time and space into account which suggest performing transient flow analysis in 3D. Furthermore, the depth and anchorage of the anti-filtration barrier are to be optimized in a parametric analysis, in which the calculations are carried out many times for different geometrical dimensions of the barrier. Thus, it is worth considering extracting smaller sub-areas from the whole area and treating them as separate numerical 3D models (computational domains). Such an approach enables to reduce time of calculations. Regardless of that, it should be noted that it may also be necessary to carry out steady-state flow analysis. It may be required due to the national regulations and standards. Such additional analysis also makes sense if a transformation of the object into a reservoir with constant water level is considered as an option for the future operation scheme.

The assessment of dam and subsoil deformation as well as the stability analyses (using the SSR method (Cała, 2007; Griffiths & Lane, 1999)) should be carried out taking into account the hydromechanical coupling. In addition, the elastic-plastic constitutive relations must be assumed for the subsoil and dam body. The Mohr-Coulomb model is most often adopted. Consequently, the stress history needs to be included in the analysis, in particular by modelling the successive stages of the dam construction. Furthermore, the factor of safety (FOS) during the flood event variates in time due to the transient groundwater flow caused by non-stationary water table level. Thus, to assess minimum value of FOS, which determines the design solutions, the stability analysis has to be performed at least a dozen times for selected moments of the analysed process. These issues entail a significant increase in computation time compared to the analyses described in step (2). Therefore, within the proposed approach, it is postulated to perform deformation and stability as 2D problems. Their domains can be extracted as the cross-sections of the three-dimensional model described in step (1). The number of cross-sections and their locations should be determined on the basis of the dam dimensions, arrangement of soil layers and terrain geometry. It is worth noticing that stability results obtained from 2D analysis are usually more conservative than 3D ones (Fredlund et al., 2010). Therefore, such an approach establishes a safe estimate.

Assessments of the deformation as well as the stability are carried out for all phases of construction and operation of the facility, in particular: after the completion of dam rise and during the process of flood detention according to the assumptions predefined in the design. In the case that the condition of overall dam stability is not fulfilled, changes in the shape of the dam, a reinforcement of the subsoil or an increase of the depth of the barrier should be considered. In all these situations, it is necessary to repeat the calculations described in steps (2) and (3).

According to the step (1) of the proposed analysis algorithm, first, it was necessary to create a 3D model of the subsoil and the dam. An area of the analyses (shown in Fig. 2) was selected based on the design documentation (Sweco Consulting Sp. z o.o. 2019).

Figure 2

As a consequence, 3D terrain surface model was created, based on the map presented above. Since the calculations were performed in the framework of FE, the next step was to import the terrain surface into the FE software. ZSoil software (Truty et al., 2020) was used in the particular case. Due to the lack of ready-to-use tools, the import was done using Python programming tools (Van Rossum & Drake Jr, 1995). Resulting 3D model of the dam substrate is presented in Fig. 3.

Figure 3

The model is further complemented with a spatial reconstruction of the subsoil layers, i.e. appropriate materials and the corresponding mechanical parameters are assigned to finite elements. The data used for the subsoil reconstruction were acquired from the boreholes profiles obtained from geotechnical documentation. Location of boreholes taken into consideration, together with the dam body, is shown in Fig. 4. The boreholes adopted into FE model are presented in axonometric view in Fig. 5.

Figure 4

Figure 5

On the basis of the reviewed documentation, it was concluded that directly under the dam quaternary formations of silty and sandy loams as well as sands and gravels with different consistency and density indices are located. Under the layer of quaternary formations, there is a rock mass built of marls, characterized by a very high variability of strength and deformation parameters. The geotechnical parameters of the quaternary soil layers were adopted according to the geological documentation. Rock mass classifications, namely Rock Quality Designation (RQD) and Rock Mass Rating (RMR), were additionally used for the estimation of the rock mass strength and deformation parameters (Bieniawski, 1988; D. Deere, 1988; D. U. Deere & Deere, 1989). A view of the subsoil layers arrangement, obtained using kriging technique, is presented in Fig. 6.

Figure 6

The analysed area with dimensions of 1200 × 400 × 100 [m] (see Fig. 6) was divided into two separate models of similar size. Such a division allows for the generation of a denser finite element mesh and results in more precise analysis of filtration phenomena within the abutments. The geometry of the model along with the reconstruction of the probable system of geotechnical layers of the subsoil is presented in Fig. 7. The results are presented only for the western abutment throughout the paper.

Figure 7

The filtration phenomena were analysed as a transient groundwater flow problem for the case of flood detention simulation as well as steady-state problem assuming the highest water level in the reservoir as a constant damming level. The boundary conditions were defined as follows: on the vertical lateral surfaces and on the bottom surface, the second-type boundary condition with zero value of the outflow through the boundary (i.e. the impermeable edge) was prescribed. The first-type boundary condition was imposed on the ground surface in the reservoir (upstream) with the value of the hydrostatic pressure corresponding to the damming level, as illustrated in Fig. 8. The remaining part of the top surface of the model is impermeable.

Figure 8

According to Section 2, the HYD limit state was verified due to the possibility of suffusion. Verification was based on the calculation of spatial distribution of the F_i parameter (Eq. (3)). An example of the distribution of the F_i parameter is presented in Fig. 9.

Figure 9

The minimum depth of anti-filtration barrier under the central part of the dam as well as its anchorage in the abutments (length and gradual reduction in depth) was established from the parametric analyses based on several conditions. First of them, given by Eq. (3), requires F_i parameter to be less than 1.0 within the entire domain at any moment of analysed groundwater flow process. Another condition was defined as the necessity to cut off the privileged filtration paths through the layers of high permeability. This requirement can be verified on the basis of spatial arrangement of the geotechnical layers in the subsoil obtained from the interpolation created using kriging technique. Permeable layers, characterized by the filtration coefficient k > 1.0e-3 m/s, as well as the geometry of anti-filtration barrier are presented in Fig. 10. Furthermore, it has to be noticed that the water discharge was calculated within the parametric analysis and the final geometry of the diaphragm (see Fig. 10) was determined so as to reasonably limit water flow under the dam.

Figure 10

The deformation analysis, leading in particular to the determination of the dam base settlement as well as the stability analysis, was carried out on 2D models. An exemplary computational model, created in the ZSoil software, is presented in Fig. 11.

Figure 11

Plane strain, full poromechanical coupling, Mohr-Coulomb failure criterion and van Genuchten model (Labuz & Zang, 2012; Van Genuchten, 1980; Xiang-Wei et al., 2010) were implemented in the FE analyses. The model consisted of 12,541 finite elements of the EAS type (Truty et al., 2020), with the denser mesh in the area of the anti-filtration barrier.

The dam erection was simulated in the model as a staged, layer-by-layer construction. The consecutive stages starting from the reference state and the initial excavation are shown in Fig. 12.

Figure 12

Displacement boundary conditions were assumed as follows: horizontal displacements have been fixed on the vertical edges of the model boundary, while horizontal and vertical displacements were fixed on the bottom edge (see Fig. 11). The edges of the model marked in purple represent so-called seepage elements. Their function is described as follows. First-type boundary condition is set if the actual pressure is prescribed. Otherwise, the action is dependent on the degree of saturation s in the seepage element: if s < 1, then the impermeable boundary is set else free outflow is allowed and the pore pressure is set to 0 on the boundary. The initial condition was defined by setting groundwater table level in accordance with geological documentation.

The flood event was simulated in a manner analogous to that used for water flow analyses (see Section 3.2). Therefore, the boundary conditions governing the water flow were assumed as the pore pressure equal to hydrostatic pressure on the bottom of the reservoir and slope of the dam. Note also that the aforementioned hydrostatic pressure has to be taken into account as a load at the model boundary in the deformation process. This time-dependent load is presented in Fig. 13 for selected moments of time during simulated flood event.

Figure 13

The calculations of the factor of safety (FOS) were performed using the SSR method (Cała, 2007; Griffiths & Lane, 1999) for selected moments of time during flood wave passage. The ZSoil software enables a stability evaluation to be carried out at any stage of the problems that are non-stationary with regard to water flow. Due to this, it is possible to determine the evolution of the factor of safety in time during the passage of a flood wave, which is related to the variable height of water damming in the reservoir. Time t = 0 h corresponds to the situation directly after completing the construction of the dam. The time t > 0 h corresponds to the simulation of the flood wave passage (transient flow). An exemplary diagram of the FOS evolution is shown in Fig. 14. This diagram also shows the FOS values with regard to the depth of the anti-filtration barrier, representing the results of the parametric analysis. Additionally, Fig. 15 shows the potential failure mechanisms determined for two time moments, in which: (a) the water level in the reservoir reaches its maximum height; (b) the lowest FOS value is obtained. Fig. 16 shows vertical and horizontal deformations of the dam and the subsoil for the case (a), specified above, i.e. for the moment corresponding to the maximum damming level.

Figure 14

Figure 15

Figure 16

The subsections above present selected results of the analyses to a limited extent. The discussion in this section summarizes the results of all the analyses performed for the facility under consideration. The results obtained are synthetically described below.

The primary purpose of the analyses carried out for the particular facility was to optimize the depth of anti-filtration barrier. It was performed based on the parametric analyses concerning: –

Hydraulic gradients

The parametric analysis consisted in the comparison of results obtained for the following pre-assumed lengths of the barrier: 0 (no barrier), 10, 13, 14, 15, 20 [m]. The outcomes of the computations showed that the values of hydraulic gradients are significantly smaller than critical values regardless of the barrier length (see Fig. 9). In other words, the critical gradient condition is satisfied in the entire domain of substrate and the dam. With regard to the overall stability of the dam, it can be concluded that the length of anti-filtration barrier does not affect the values of FOS except of the range of time (20–70 hours) corresponding to high water level, where minor fluctuations are observed (see Fig. 14). In order to minimize the discharge of water under the dam, 3D spatial arrangement of subsoil layers, obtained using kriging technique, was examined. It was assumed that the barrier must cut off the layers of high permeability, i.e. characterized by the filtration coefficient k > 1.0e-3 m/s (see Fig. 10). This requirement determined the final depth of anti-filtration barrier (as equal to 14 m) as well as the length of its anchorage in the abutments.

Stability and deformation analyses were conducted for four dam cross-sections in the case of a dry reservoir (after dam construction completion) as well as during a flood event. Maximum settlement of the subsoil in the level of dam foundation (resulting from dam construction) occurs below the crest. The maximum settlement value was evaluated as 11.3 cm. Vertical displacements caused by the flood event can be described as downward movement in the area of upstream slope and uplift (upward movement) in the area of dam crest and downstream slope. These displacements generally do not exceed 2–3 cm (see Fig. 16).

The values of FOS were determined at selected time steps in order to identify the moment corresponding to its minimum value, namely FOS_min. Fig. 14 shows an example of the evolution of FOS in time during flood event for a selected cross-section. The minimum value of FOS was identified for time t=96 h. Finally, it was stated that the stability of the dam and the subsoil in all cross-sections is preserved. For a dry reservoir, the margin of safety is not less than 1.40, while during flood event, the minimum value of margin of safety is not less than 1.34.