The article proposes the methodology of designing dams of dry flood control reservoirs. The algorithm is developed so as to meet all the requirements given in the Eurocode 7 and, at the same time, to be efficient in terms of necessary calculation time. Furthermore, the presented numerical procedure enables the optimization of design solutions, e.g. the depth and length of the anti-filtration barrier, by means of parametric analyses. The approach assumes the use of numerical methods, in particular, finite element (FE) analysis. Three-dimensional (3D) reconstruction of the terrain topography and subsoil layer arrangement performed in step (1) sets the base for further analyses. In step (2), the filtration phenomena are assessed based on the 3D analysis of a transient groundwater flow. In step (3), the state of displacement is evaluated and the stability is verified for all the relevant phases of construction and operation of the facility, in particular, in the course of simulated flood detention. The analyses in step (3) are carried out on 2D models corresponding to the design cross-sections of the dam. This significantly reduces the computation time (compared to 3D analysis) and, at the same time, provides a safe estimate of factor of safety. The performance of the proposed algorithm is shown on the numerical examples of the computations concerning the dam of Szalejów Górny dry anti-flood reservoir located in Poland.

#### Keywords

- finite element method
- stability
- settlement
- filtration
- kriging

Sudden, increased inflow of water to the river may cause water to overflow the riverbed and flood adjacent areas causing significant damage of urban as well as agricultural areas. In such events, financial losses often exceed billions of dollars as in 1997 when during the floods in Poland and the Czech Republic damages were estimated at $3 billion (Kundzewicz et al., 1999; Szamalek, 2000). In order to avoid or reduce the amount of damages caused by flood, it is necessary to build dedicated anti-flood infrastructure. Therefore, dry anti-flood reservoirs present great potential in controlling the floods. The main role of dry anti-flood reservoirs is short-term retention of floodwaters for a time needed to maintain a safe downstream water flow (Matyas et al., 2009; Poulard et al., 2008). An important feature of this type of construction is that the reservoirs are empty most of the time and retain water only in the case of flood risk (Łydżba et al., 2017). Compared to classical multipurpose reservoirs, which are filled with water for a lifetime, the construction of dry reservoirs is definitely cheaper, and they also have less impact on the environment. Beyond the flood event, the reservoir area is usually used as a grassland or meadow. If there is no flood risk at a given time, the river flows freely through the concrete conduit of the dam body. If a flood occurs, water is dammed up in the reservoir. The outflow from the reservoir is then being controlled so that to keep the downstream areas safe. Eventually, in an extreme situation of reservoir overfill, a spillway starts to convey excess of the water so that to ensure the safety of the dam structure itself, mainly the stability of its body. Dam of dry anti-flood reservoir can be either made from soil or if needed of concrete (Matyas et al., 2009). The choice of the material to be used for construction depends on dam location on the river, the width of the valley and the availability and ability to deliver the material to the construction site. Earthen dam usually has a core made of cohesive soil within the body. Under the core, an anti-filtration barrier is installed. The main purpose of the barrier is to extend the filtration path and, thus, to limit the underground water flow below the dam. Regardless of the dam type, such a structure must be equipped with the necessary conduit. It is used to safely conduct the water flowing downstream through the dam. In case of dry anti-flood reservoir, the conduits are generally unenclosed and have a width similar to the riverbed. The dimensions of the reservoir dam, i.e. its height as well as the number of conduit devices, their size and location should be determined based on hydrological analyses.

Properly designed reservoirs are an important part of flood control system. History has shown that dry anti-flood reservoirs have saved large areas from damage, e.g. in the United States, major floods were avoided on the Miami River in 1937, 1959, 1982 and 1991 (Matyas et al., 2009). On the other hand, ensuring the safety of a dam structures itself emerges an extremely important issue due to the dynamic nature of possible failures, which can result in large numbers of fatalities and significant financial losses. The Banqiao Dam failure in China, for example, resulted in the deaths of more than 26,000 (Xu et al., 2012), or losses from the Teton Dam failure in the United States have been estimated at $2 billion (Bolton Seed & Duncan, 1987; Panel, 1976). Considering the possible risks, the design of a dry anti-flood reservoir dam should be based on a comprehensive approach that analyses all potential failure modes.

In the European Union (UE), the applicable standard for the design of geotechnical structures such as earthen dry anti-flood reservoir dams is Eurocode 7 (

For the development of the design procedure, the requirements of Eurocode 7 were adopted as the prerequisite. Therefore, the following ultimate limit states should be computationally verified regarding the earthen dams: failure or excessive deformation of the soil substrate due to loss of overall stability (GEO), internal failure or excessive deformation of the structure, e.g. conduit tunnels or any reinforcing elements such as retaining walls and anchors, if present (STR), effects of the buoyant force of water (UPL) or unfavourable phenomena related to water flow leading to hydraulic uplift of soil particles (HYD). Moreover, the Eurocode 7 states that the analysis of hydraulic gradients, pore water pressures and filtration forces should take into account the variability of soil permeability over time and space, the spatial variability of water levels and pore water pressure over time and any changes in boundary conditions.

The aforementioned requirements follows, in particular, that the analysis should be performed under the assumption of transient flow. Furthermore, to take the variability of soil permeability over space into account in an unequivocal manner, the three-dimensional models should be used. This states a complex numerical problem to be solved. Consequently, the algorithm proposed in this work was developed in the framework of FE analysis and consists of the following 3 general steps.

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:
_{cr}_{i}

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

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:
_{i}_{i}_{i}_{i}

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).

The method presented in previous section was used to carry out verification calculations for the dam of the Szalejów Górny dry anti-flood reservoir. The facility is located on the Bystrzyca Dusznicka River in south-western Poland. The reservoir is constructed within the framework of the Odra-Vistula Flood Management Project which is co-financed by the World Bank funds (The World Bank, 2021). The primary aim of this programme is to increase the safety level of anti-flood protection in the large areas of Poland. Within this project, new flood control reservoirs are to be built in the Kłodzko Valley, inter alia: Boboszów, Roztoki Bystrzyckie, Krosnowice and Szalejów Górny.

The considered Szalejów Górny reservoir has the maximum capacity of 9.9 million m^{3} and the flooding area at maximum water level of 118.7 ha. The reservoir controls 64% of the entire Bystrzyca Dusznicka River basin, which significantly affects anti-flood protection of the downstream areas, including the town of Kłodzko. The length of the dam in the crest axis is 735 m, and its height is up to 19.3 m. A characteristic cross-section of the dam is shown in Fig. 1.

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).

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.

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.

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.

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.

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.

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 _{i} parameter (Eq. (3)). An example of the distribution of the _{i} parameter is presented in Fig. 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 _{i}

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.

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.

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

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.

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

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

Overall stability of the dam

Average discharge of water under the dam

Spatial arrangement of geotechnical layers under the dam

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

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.

The paper proposes a comprehensive approach to the design of earthen dams of dry flood protection reservoirs. For the development of the design procedure, the requirements of Eurocode 7 were adopted as the prerequisites. As a consequence, all the required limit state conditions are verified so as to meet all standard requirements and at the same time to reduce the calculation time to the necessary minimum. It is, furthermore, postulated to optimize the design solutions such as depth and length of the anti-filtration barrier, dam slope steepness or strengthening of subsoil by means of a parametric analysis. Detailed parameters and types of analyses can be adapted to the specific conditions of the objects being designed. Consequently, the design optimization enables to reduce construction costs while maintaining an appropriate safety margin.

The approach assumes the use of finite element method. Three-dimensional reconstruction of the terrain topography and subsoil layer arrangement performed in step (1) which sets the base for further analyses. In step (2), the filtration phenomena are assessed based on the 3D analysis of a transient groundwater flow. At this stage, the HYD ultimate limit state is verified with regard to the possibility of suffosion. Additionally, at this point, the required geometry of the anti-filtration barrier is optimized. For this purpose beside the condition of HYD limit state condition, the flow rate under the barrier (and next to it in the area of abutments) is taken into account. Furthermore, it is postulated to cut off the filtration paths through the permeable layers of the substrate with the barrier. The verification of general stability analysis (GEO) as well as the filtration stability (HYD and UPL) is performed in step (3). At this stage, the effects of hydromechanical coupling are taken into account, which generally extends the calculation time. Therefore, it is postulated to perform calculations in step (3) for selected computational cross-sections extracted from the previously created 3D reconstruction of the substrate. The number of sections and spacing between them (along the dam axis) depends on dam dimensions and the variability of geotechnical conditions in the substrate. The stability as well as the state of displacement is evaluated for all the relevant phases of construction and operation of the facility, in particular, in the course of simulated flood detention.

The proposed method was successfully used for the design purposes of several dry flood control reservoirs in the Odra basin in the south-western part of Poland. The article presents the selected results of numerical calculations performed in ZSoil software concerning the dam of the Szalejów Górny dry flood control reservoir.

Seismic bearing capacity of shallow strip footing embedded in slope resting on two-layered soil FEM modelling of the static behaviour of reinforced concrete beams considering the nonlinear behaviour of the concrete Evaluation of Tunnel Contour Quality Index on the Basis of Terrestrial Laser Scanning Data Efficiency of ventilated facades in terms of airflow in the air gap Characteristic parameters of soil failure criteria for plane strain conditions – experimental and semi-theoretical study A comprehensive approach to the optimization of design solutions for dry anti-flood reservoir dams Laboratory tests and analysis of CIPP epoxy resin internal liners used in pipelines – part II: comparative analysis with the use of the FEM and engineering algorithms Usefulness of the CPTU method in evaluating shear modulus G changes in the subsoil_{0}