A comprehensive computer aided approach to design of earthen dams of dry flood-control reservoirs

A generalized, comprehensive approach to geotechnical design of earthen dams of dry flood-control reservoirs is developed. It consists of three steps:• In Step 1 three-dimensional (3D) model of terrain and dam subsoil is created. The arrangement of geotechnical layers is reconstructed based on results of geotechnical investigation using geostatistical interpolation, namely, kriging.• Step 2 involves 3D finite element (FE) analysis of transient groundwater flow through the subsoil. Anti-filtration barrier extent is determined by appropriate parametric analysis and the condition of critical hydraulic gradients is verified.• Step 3 concerns deformation assessment and stability analysis. Computations are carried out for selected cross-sections of the dam. At this stage transient coupled problem of deformation and groundwater flow is considered. Stability is assessed with the use of shear strength reduction technique. The proposed approach enables verification of all relevant conditions of limit states. Furthermore, it has been developed in such a way as to enable computationally efficient analysis of considered phenomena. The method constitutes a general framework of design procedure. It can be easily adapted to conform any standard requirements, e.g., by taking adequate values of partial safety factors. Moreover, design solutions can be optimized through parametric analyzes.


Method details
Overview Dry flood-control reservoirs are important elements of passive flood protection infrastructure. Such reservoirs remain empty most of the time and fill up during high water [1] . This way, the outflow from the reservoir is limited so that to keep the downstream areas safe. To ensure proper and safe operation of such facilities, they must be properly designed. The design method presented in the article concerns geotechnical aspects and applies in particular to earthen dams of flood-control reservoirs. From technical point of view the most important issue is to ensure the stability of dam and its substrate with regard to mechanical strength as well as phenomena related to groundwater flow. A significant difficulty in designing this type of structures is the fact that they temporarily dam up water. This entails transient states of groundwater flow. As a further consequence mechanical stability (in terms of factor of safety) changes over time in the course of flood detention. Moreover, the water filtration in the substrate is intrinsically three dimensional (3D) and this should be taken into account in the design procedure. Thus, the proposed method is based on creation of 3D models followed by appropriate numerical computations in order to verify the stability of the dam and conditions of safe water flow in the substrate. Originally, the method [2] has been developed so that to fulfill requirements given in Eurocode 7. Nevertheless, it can be easily adapted to conform any standard mostly by taking into account the appropriate safety factors. The presented method has been formulated in such a way as to allow credible analysis of the phenomena occurring within the dam and its substrate. It assumes using numerical methods, most preferably finite element (FE) analysis, to obtain solutions of initial-boundary value problems. They set the basis for verification of appropriate limit states. The procedure assumes rational way of using numerical models with minimal use of resources to reduce the computation time. Simplifications of the model are proposed in order to reduce the 3D computational domain of groundwater flow analysis. Mechanical stability is assessed using shear strength reduction (SSR) technique. Since it requires application of hydromechanical coupling which introduces significant complexity to the model the task is implemented for extracted cross-sections as 2D problems. Furthermore, the proposed framework of computation not only enables verification of appropriate conditions of limit states, but also adjustment of design assumptions, e.g., an extent of the cut-off screen. The method is illustrated with examples of its application for two designed recently dry flood-control reservoirs located in south-western Poland, i.e. Szalejów Górny [3] and Krosnowice [4] . Computations for the mentioned objects were carried out using ZSoil finite element method code. However, the use of proposed method is not restricted to a particular type of software.

Step 1. Reconstruction of the terrain surface and probable soil layers arrangement
Designing a dam for an earthen flood-control reservoir requires a sequence of computations performed on a three-dimensional numerical model. However, the creation of such a model must be preceded by specifying the area covered by the analysis, mapping the actual surface of the terrain and reconstructing the probable arrangement of geotechnical layers. These three tasks are carried out on the basis of the available design documentation containing a hypsometric map of the terrain and results of geotechnical investigation, e.g. boreholes or soundings. The domain should cover at least the area of the dam. However, to minimize the influence of boundary conditions on the results of numerical calculations, it is strongly recommended to extend the domain by at least the total width of the dam in each direction ( Fig. 1 ).
A three-dimensional terrain mesh model can be generated using e.g. CAD tools. It is then transferred into the FE software (where the numerical simulations are performed) and creates the top boundary of the 3D computational domain. The model must also be bounded from below, and it is important here that the depth of the subsoil should be at least several times greater than the height of the dam. The exemplary preliminary models of the terrain and subsoil, together with the division into finite elements are shown in Fig. 2 .
The final numerical model is created by supplementing the preliminary one ( Fig. 2 ) with the reconstructed system of geotechnical layers. In practice, the reconstruction consists in assigning the appropriate material and its parameters to each element. This task can be performed with the use of various techniques, however, geostatistical interpolation, e.g. kriging [5] , seems to be the most rational here. The interpolation treats the elevation of layer bottom and ceiling as a random variable [ 6 , 7 ]. Its values are fixed in the locations of available boreholes and interpolated in the entire domain. The number of available boreholes depends on the implemented subsoil investigation program. It is obvious that the more information about the subsoil is available, the better the reconstruction of the actual layer arrangement under the dam. Usually, for economic reasons, the number of boreholes or soundings is limited to the minimum requirements in accordance with applicable standards and regulations. In Fig. 3 available boreholes (treated as fixed data within kriging) for considered dams are presented.
Interpolation by kriging requires the assumption of the so-called semi-variogram function which measures "the strength " of statistical correlation as a function of distance. In this regard, the spherical and exponential models are most often used [8] . An important parameter to be established is also the correlation length/radius [6] . This can be evaluated by fitting models to sample semi-variograms by weighted least squares approximation or by assuming value of correlation length in relation to the dimensions of the dam [ 8 , 9 ]. Utilizing the given data provided by available boreholes ( Fig. 3 ) the process of interpolation consists in minimizing the error variance of the interpolated variable. Fig. 4 presents numerical models with reconstruction of soil layers spatial arrangement obtained with spherical models and assumption that correlation length is equal to 500 m. The result of kriging is shown in part of the volume in the back, while in the foreground only the location of boreholes is presented.

Step 2. Evaluation of groundwater flow and determination of cut-off screen depth
Creation of the numerical model as described in Step 1 is followed by groundwater flow assessment. In particular, cut-off screen depth and its anchorage length in the abutments are determined by a proper parametric analysis. The extent of anti-filtration barrier must be enough to meet the conditions described next. The main consideration in this context is the limit state related to unfavorable phenomena resulting from groundwater flow in the subsoil such as internal erosion, soil boiling or soil uplift. The conditions of soil boiling or uplift occurrence are verified during the numerical analysis of the overall dam stability (see Step 3). To avoid possible occurrence of internal erosion, i.e. suffosion, the critical filtration gradients must not be exceeded. It is expressed by the following formula: where: i -hydraulic gradient, i cr -critical values of hydraulic gradient for a given soil, i -safety factor determined according to national regulations and/or relevant standards. In order to verify condition (1), first an initial-boundary value problem concerning transient flow in partially saturated medium has to be solved. In this regard, commercial software most often utilize van Genuchten model of soil-water retention curve [10] . Furthermore, Richards hypothesis of gas phase continuity in the aeration zone allows one to perform single phase flow analysis. Such an approach is recommended in the proposed method. The computations are performed using the 3D model resulting from Step 1. For saving computational resources the entire model can be divided into two parts. Then, each of two abutments are considered   separately. This routine is optional. However, it was utilized in the examples illustrating the method. Appropriate models are depicted in Fig. 5 . For a given model domain boundary conditions (BCs) are to be defined. The groundwater table level known from geotechnical investigation is extrapolated to the vertical planes limiting the model domain on upstream and downstream sides. The relevant BCs should be set as pore water pressure linearly increasing with depth. At the level of water table the value of pore pressure is set to 0 and a constant gradient of the pore pressure along elevation ordinate corresponds to the hydrostatic increase. Remaining parts of model boundary are at first set to be impermeable. For boundary conditions set in this manner the steady state problem of groundwater flow is solved. Obtained distribution of pore pressure is set as initial condition for the following part of computation involving analysis of transient filtration flow. In order to perform such analysis a new boundary condition is to be defined on the part of a boundary which corresponds to the basin of the reservoir. This BC represents hydrostatic water pressure exerted on the basin in the course of flood wave attenuation. This BC is depicted in Fig. 6 . The distribution of actual water pressure set as the boundary condition results from the anticipated, time-varying damming level in the reservoir. In general the function of water table level in the reservoir versus time is calculated based on the basin geometry and stream discharge function assumed for the design.
The distribution of the hydraulic gradient i is determined according to Darcy's law based on the fluid velocity distributions obtained from FE analysis. Within the saturation zone the following formula applies:  where: − hydraulic conductivity of soil, − fluid velocity.
The critical values of hydraulic gradients i cr in different soil types are determined most often from empirical relationships. Among the most widespread there are the formulas of Sichardt [11] and Abramov [12] , with the latter being more conservative. Eventually, the Abramov formula is proposed to be applied. In terms of critical gradient it expresses as: where k is given in m/s. It is inconvenient to verify condition (1) using its direct form because the critical values of the hydraulic gradient are generally different for each layer. For this reason, it is proposed to transform the condition (1) to the following form [1] : The parameter F i is the multiplicity of the limit admissible gradient. Then, the interpretation of It is to be noted that in general the condition (3) locally can be violated without indicating dangerous water flow states. This situation may take place at the base of the cut-off wall due to the singularity of numerical solution next to the "sharp edge" of subdomain. In addition, minor, local violations of condition (3) may be conditionally accepted as long as they are below the ground surface, or next to it with flow velocity vector directed to the ground.
Furthermore, it is advised to place the cut-off wall so that to partition permeable layers to limit the filtration flow under the dam as well as through abutments. It is convenient to visualize the permeable layers in 3D model to inspect whether they are split with the anti-filtration barrier. Such views of the model are presented in Fig. 8 . The particular layer is considered permeable if its hydraulic conductivity k > 1.0e-3 m/s (see blue layers in Fig. 8 ) or it is at least two orders of magnitude larger than hydraulic conductivity of  a dominant layer in the subsoil at the level of cut-off wall base. Such a case is marked orange in Fig. 8b . If permeable layers cover larger areas of the subsoil then the extent of cut-off wall should be determined in parametric analysis by testing different depths and anchorage lengths and checking whether the condition (3) is met. In this context, please note considered variants of cut-off wall anchorage lengths marked with different colors in Fig. 8 .
Above-described requisites regarding groundwater flow states are referred to as necessary conditions. These must be considered in conjunction with the further requirements described in the next step. Moreover, it may also be necessary to carry out steady-state flow analysis for maximum damming level assumed in the design. It may be required due to the national regulations or 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.

Step 3. Deformation assessment and stability analysis
Analyses within the scope of Step 3 for the assessment of deformation and settlement of subsoil, as well as stability analyses (using the SSR method) [ 13 , 14 ] are carried out for selected plane cross-sections of the dam. 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 location is to be determined according to the size of the dam, the arrangement of the soil layers and the geometry of the terrain. It should be noted that the stability results obtained from a 2D analysis are generally more conservative than from a 3D analysis [15] . Therefore such an approach establishes a safe estimate. The elastic-plastic constitutive relations must be assumed for the subsoil and dam body. The Mohr-Coulomb model is adopted in this regard. Furthermore staged construction of the dam is taken into account. In other words, erection of the dam is modeled as a stepwise (sequential) "extension " of the model geometry by successively "appearing" layers. Figs 9 and 10 show examples of FE numerical models of dam cross-sections and selected construction stages included in the analyses.
Settlements resulting from the construction of dam are computed as vertical displacements at the foundation level. Unless it is automated in the particular software used the settlement is to be calculated as the increment of vertical displacement between reference state (before dam construction) and the one corresponding to the completion of the structure. Examples of calculated settlement are presented in Fig. 11 . With regard to the settlements calculated in this way, the condition of the appropriate limit state is met if the maximum value of obtained displacement is not greater than maximum admissible value. The latter is given in the relevant standard or legal requirements applicable in a given country.
The computations of dam subsoil deformation as well as stability assessment is carried out taking into account the hydromechanical coupling. This entails a significant increase in computation time compared to uncoupled analysis described in Step 2. This is the main reason for performing deformation and stability analyses as two-dimensional problems within the proposed design approach. On the other hand, using so called effective parameters of strength and considering coupled problem enables to verify a few limit states at ones. It includes general stability as well as filtration related phenomena that may occur, i.e. soil boiling and soil uplift mentioned earlier. All these limit states require the same condition to be met. The factor of safety ( FOS ) has to be not less than FOS min . The latter is minimum required value of FOS provided by relevant standard or legal requirement. Regardless of particular value of FOS min (and any partial safety factors) it can be convenient to represent condition of stability in the form: where F ( t ) stands for relative safety margin. Due to the transient groundwater flow F ( t ) varies in time in the course of flood event. Thus, in order to evaluate the minimum value of FOS on which the design solutions depend, the stability analysis should be carried out at least a dozen times for selected moments in time of the analyzed process. The examples of F ( t ) changes with time are presented in Fig. 12 . Note that for different objects considered variation in time is either apparently pronounced (as for Szalejów Górny reservoir -see Fig. 12a ) or not present (see Fig. 12b related to Krosnowice reservoir). Fig. 13 shows the potential mechanisms of failure for different moments of time. Assessments of the deformation as well as the stability are carried out for all the construction and operational phases, in particular after completion of the dam construction and during the damming process according to the assumptions made in the project. In the event of non-compliance with the limit states, changing the shape of the dam, a reinforcement of the subsoil or an increase of the depth of the barrier have to be considered. In such situations it is necessary to repeat the calculations described in Steps 2 and 3.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The data that has been used is confidential.
involved in works carried out on the basis of commercial orders from Hydroprojekt Wroc ł aw Sp. z o.o. and Sweco Consulting Sp. z o.o. in the period of 2014-2020, e.g. [ 3 , 4 ]. The method reported in the present article has been formulated using the experience gained during the execution of the above-mentioned works.
This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.