Reliability analysis of an earth dam in operating conditions using direct coupling

In this study, we present a numerical investigation of the direct coupling between the deterministic GeoStudio package (Seepage/W and Slope/W software) and the StRAnD reliability package for normal operating conditions of an existing old earth dam over time. Direct coupling avoids errors associated to point estimate methods and to using response surface surrogates. One relevant feature of the study is the use of realistic pore water pressures for each equilibrium analysis, accomplished by considering a long-term steady-state analysis as an initial condition, followed by four years of equilibrium analysis, before each transient seepage analysis. All reliability analysis were performed using the first-order reliability method. The most important random parameters found in sensitivity analyses are four dam body parameters (saturated hydraulic conductivity, ks; specific weight, γ; effective cohesion, c′; and friction angle, ϕ′) and two filter parameters (ks and ϕ′). A range of values of the relationship between the reliability index (β) and the factor of safety (FS) was found for all probabilistic and deterministic results. Finally, a large difference in terms of critical deterministic and probabilistic slip surfaces is identified for the same time of analysis. Realistic pore water pressures used in dam equilibrium analysis, by considering random seepage analysis in previous 5 to 30 days. Direct coupling of deterministic and reliability softwares avoids errors associated to point estimate and response surface surrogate methods. Existing old earth dam studied in Normal Operating Condition. For same analysis time, large differences are observed between critical deterministic and probabilistic slip surfaces. Realistic pore water pressures used in dam equilibrium analysis, by considering random seepage analysis in previous 5 to 30 days. Direct coupling of deterministic and reliability softwares avoids errors associated to point estimate and response surface surrogate methods. Existing old earth dam studied in Normal Operating Condition. For same analysis time, large differences are observed between critical deterministic and probabilistic slip surfaces.


Introduction
A country's economic growth is directly related to the construction of dams. The intended purposes of large dams usually include providing water for irrigation, water supply to cities, improving navigation, generating hydroelectric power and flood control [1]. Few dams [2] serve all these purposes (e.g., Peixotos, Taiaçupeba, Capivari-Cachoeira, Euclides da Cunha Hydroelectric Power Plants in Brazil), but some multi-purpose dams serve more than one (e.g., Itaipu, Belo Monte, Tucuruí, and Santo Antônio Hydroelectric Power Plant in Brazil).
The failure of large dams is a concern in many countries due to the high economic and social consequences associated with it. Zhang et al. [3] present a study of failures in dams from over 50 countries; additionally, Menescal [4] lists 166 cases of accidents and incidents that occurred in Brazil from 1954 to 2009. In a study of a new dam or during the construction, the geotechnical engineer must be able to make reliable predictions about the behavior of the dam, under every anticipated operating condition (e.g., end of construction, first filling of the reservoir, rapid drawdown and normal operating conditions). The predictions usually involve both judgments and quantitative analyses, based on data provided by the site investigation team [2]. All the information about the dam must be available (e.g., topography, hydrology, geology, laboratory and field test, …, etc.), and additional laboratory or field tests (e.g., geophysics, additional perforations, …, etc.) may be required to perform a complete study.
The safety evaluation of aging dams is an essential but complex undertaking, especially when the original foundation investigation, dam design and construction details are not known or are associated with significant uncertainty [3,5]. When information is limited, the real condition of the dam can often only be estimated. Continuous monitoring of the dam by the variation in the water reservoir is an important activity of dam safety management during normal operating conditions (NOC).
In a deterministic approach, the safety condition of a dam is usually verified by Factors of Safety (FS). During the last four decades, different authors [6,7] have worked to quantify the variability in geotechnical materials and in developing methodologies to perform reliability analysis of slopes and embankments in different conditions [8][9][10][11]. Kumar et al. [12] performed reliability analysis in a zoned dam using a deterministic GeoStudio package (Seepage/W and Slope/W software) with two different methods specifically multivariate adaptive regression splines (MARS) and relevance vector machine (RVM). Reliability analysis of earth dams using response surface methodology, in combination with first order reliability method and numerical analysis was presented by Sivakumar Badu e Srivastava [13]. Additionally, the random field theory is emerged in the field of geotechnical during the last years in contrast to the random variables. Guo et al. [14] performed different types of random fields in a real earth dam. An advanced machine learning algorithm called XGBoost to evaluate the earth dam slope failure probability in spatially variable soils was developed by Wang [15]. Other different studies [16][17][18] use seepage and stability analysis in dams to perform stochastic analysis with random finite elements. These approaches are different in mechanical terms and it takes more computational effort.
In new dams or dams with less than 5 years of operation, the dam is analyzed from empty to the actual state of the dam. This is possible using available information, and the computational cost is admissible. However, in old dams or dams with more than 5 years of operation (in our case, more than 50 years of operation), the dam is analyzed with a different approach. If monitoring information is available, the best approach is performing a long-term steady-state analysis as the initial condition, followed by some years of equilibrium analysis (in our case, 4 years as a recommendation). The numerical pore water pressures are calibrated based on dam instrumentation readings. The numerical calibration gives reliable results.
In this paper, an existing old earth dam is used as study object in a reliability analysis of NOC. A deterministic software (GeoStudio 2019, [19,20] and a structural reliability software (StRAnD [21] are combined to perform reliability analysis of dams. This approach is known as Direct Coupling (DC) in the literature [22,23].
A relevant feature of the study is the use of realistic pore water pressures for each equilibrium analysis, which is accomplished by considering a long-term steady-state analysis as the initial condition, followed by four years of equilibrium analysis, before each transient seepage analysis.
Recent dam failures in Brazil forced the government to change the laws and make them stricter. Different laws were passed in the last 10 years to regulate different types of dams [24]. These laws require dam owners to perform periodic evaluations of dam stability to check the current general safety state of the dam. The aim here is to provide a tool for monitoring old earth dams, considering deterministic and probabilistic approaches.
Dam safety management using a probabilistic approach represents an enhancement over traditional dam safety practices, through an integrated process that involves quantifying the uncertainties in geotechnical and operational parameters, and yielding a quantitative measure of dam safety: the probability of failure, as a function of time. An existing old earth dam in NOC is the object of study herein.
The remainder of this paper is organized as follows. The problem setting is presented in Sect. 2, which also briefly describes the performance function and reliability analysis techniques employed herein. The complete direct coupling technique for transient analysis and the characteristic of the studied dam are described in Sect. 3. The results of the transient numerical dam equilibrium analysis are presented and discussed in Sect. 4. The concluding remarks are presented in Sect. 5.

Numerical soil conditions
The Richards equation [25] represents the movement of water in unsaturated soils. This equation was obtained by the modification of Darcy′s law. The computational software SEEP/W [19] uses this formulation for solving a transient and two-dimensional seepage analysis. In the case of anisotropy, the anisotropy coefficient r k = k x /k z is defined in the function of the horizontal k x and vertical k z permeabilities. In the hypothesis of completely saturated soil, the permeability at saturation is assumed to be constant [19,26].
The Soil Water Retention Curve (SWRC) is fundamental for analyzing the behavior of unsaturated soil in a variety of geotechnical applications [27,28]. The empirical equation proposed by Van Genuchten [29] was used.
The Mohr-Coulomb criterion was extended by Fredlund et al. [30] which is available in the software Slope/W. Two independent stress state variables, namely net normal stress and soil suction are commonly used to represent the shear strength of unsaturated soils.

Random seepage analysis and random transient pore water pressures
In transient problems, the equilibrium condition changes as a function of time, mainly due to the behavior of the water reservoir (variation of the water reservoir level). These conditions produce changes in pore water pressures (PWP) which affect the equilibrium condition. In a deterministic study, the safety factor (FS) is understood to change over time, following the variation in the water reservoir level. In a similar way, in a probabilistic study, the probability of failure (P f ) and the reliability index (β) describe the safety of the dam over time. However, the point in time is unknown when the minimum FS, P f and β occurs. Thus, the deterministic and probabilistic analyses have to be performed over time using a discretized time step (t k , k = 0,1,…,K).
The reliability analysis for limit equilibrium at any time t k depends on the pore water pressures at that time. However, because some seepage parameters are random, the PWPs at time t k are random as well. Describing randomness in PWPs is difficult, because of their distribution in space. Alternatively, the correct random PWPs at time t k can be found by running the random seepage analyses since time zero [31]. This also means running the full random seepage analysis again, from time zero, for the PWPs at time t k+1 . This has a significant computational cost, when numerical solutions with many degrees of freedom are computed.
In the Normal Operation Conditions (NOC) studied herein, the reservoir water level changes constantly, going up and down, attenuating large deviations in PWPs, which could result from extreme realizations of volumetric water content or hydraulic conductivity. A practical consequence is that analyzing the whole seepage time-history is not necessary. This is very relevant in the context of this study, as the finite element models employed are very large, and the life of the studied dam is very long.
In this paper, the reliability analyses for limit equilibrium were performed considering different time intervals for the random seepage analysis (sp = 1, 5, 10 or 30 days). Hence, the reliability analysis at time t k starts at time t k-sp . In all cases, the initial condition of the dam is the mean PWP found in the deterministic (mean value) analysis. For sp = 5 days or more, no differences were observed in the calculated reliability indexes. Hence, sp = 5 days was employed in all the reliability analyses discussed herein. It means that random seepage analysis will be begun 5 days before to the desired day.

Performance function
To perform structural reliability analysis, it is convenient to describe failure events in terms of performance functions g(X), where X = {X 1 , X 2 ,…, X n } denotes the vector of random input parameters. The limit state function, g(x) = 0, sets the boundary between failure (x|g(x) < 0) and non-failure (x|g(x) > 0) states. In slope stability problems, the limit state function is expressed by the following equation [6]: where FS is the factor of safety with respect to stability, a dependent variable [32]. The critical surface and FS are calculated for each realization x during the search for the design point. The probability of failure (P f ) is given by: (1) g( ) = g X 1 , X 2 , … , X n = FS X 1 , X 2 , … , X n − 1. where f X (x) represents the joint probability density function of the random variable vector X and the integral is carried out over the failure domain. Analytical solution of Eq.
(2) only exists in particular cases of limited practical interest. In most practical applications, approximate methods must be employed, as described in reference works [33][34][35].

The first-order reliability method
The First-order Reliability Method (FORM) linearizes g(X) at the so-called design point y * in the transformed standard Gaussian Y-space, by means of a transformation y = T(x).
The y * point is found by solving the following constrained optimization problem: Find y* which minimizes y T y, subjected to g(y) = 0, where g(y) is the limit-state function in the Y-space. The smallest distance between y * and the origin of the standard Gaussian space is the reliability index:β = √ ( * ) t ( * ) . In FORM, the failure probability is obtained by approximating the limit state by a hyperplane, centered at the design point: where Φ(.) is the standard Gaussian cumulative distribution function [33].
The Hasofer-Lind Rackwitz-Fiesler (HLRF) algorithm is a very popular algorithm for solving the constrained optimization problem. The HLRF algorithm is employed in this study, but no convergence problems were encountered [33][34][35]. The gradient vector of the limit-state function is evaluated numerically via progressive finite differences method.
In this study, a coupling between the deterministic and the reliability software was developed. Each time the reliability software requires the evaluation of the limit state function, the deterministic software is called with the current set of input parameters. Then the limit state function is evaluated. This technique is called direct (DC) in the literature.
The vector that indicates the direction of the design point y * , called vector of direction cosines, is given by: where ∇g is the gradient of the limit state function with respect to the random variables. From Eq. (4), sensitivity factors 2 i are computed: these factors reveal the relative contribution of each random variable to calculated failure probabilities [33][34][35].

Transient analysis using FORM
In transient analysis, the importance of the random variables (factors 2 i ) changes over time. The FORM solution depends on the cumulative effect of uncertain seepage variables up to analysis time t k . As mentioned before, random seepage analyses were performed for sp = 5 days, starting at time t k−5 .
In contrast to FORM, in the Monte Carlo Simulation (MCS) each simulation is computed from the initial time up to final time t K , and the final solution show all the history of PWPs. The probabilistic limit equilibrium analysis at time t k can be performed by sampling stability variables, and by sampling one of the PWPs curves generated in the random seepage analysis at time t k . This is an advantage over FORM. However, the FORM analysis is still competitive if the preceding seepage analysis is limited to five days. FORM would lose competitiveness if the seepage analyses had to be performed from time zero.

Review of target reliability index standards
The target reliability indices (β target ) for different civil structures (building, bridges, etc.) are defined in various structural design codes, but this is not the case for dams. As an example, evaluation of dam safety in Brazil does not use target safety index by two major arguments: (i) the safety format used is based on safety factors; (ii) authorities do not specify the required or minimum safety levels.
Reliability indices provide a qualitative estimation of a dam performance taking into account uncertainty in loads and materials. Dams are classified by the expected performance level [36]. The target reliability values shown in Table 1 should be used in general. The β target for structural design is determined by calibrating to existing practice [37]. The target safety values of different structures cannot be used for dams. The reasons are as follows: (i) the reliability index is nominal and comparison between other structures (especially by the difference in failure modes and loading conditions) may not be correct (e.g., dams, slopes or foundations with buildings, bridges or power towers),(ii) consequences might not be comparative; (iii) Generally, dam safety risk management have large scale consequences, and this affects societal judgment.
A dam is usually designed abiding to national or international norms. Generally, a numerical model is used with information available from laboratory or field tests in critical and normal conditions of load. Computational results prove whether the modeled structure meets the safety requirements. The code requirement is: For existing dams, the following requirement should be met: To validate these two statements, different actions may be applied in the existing dam. Ditlevsen and Madsen [38] define three different categories of actions: (i) let the structure without changes,(ii) strengthen the structure; (iii) demolish the structure.

Relationship between reliability index (β) and factor of safety (FS) for specific dams
Usual safety analysis, monitoring and inspection of dams is based on safety factors. Structural reliability analysis provides a more comprehensive quantitative estimate of dam safety, but it is more complex to evaluate. Hence, it is helpful to have a measure of the relationship between reliability index and safety factors, even if this measure is valid only for very specific conditions. If such a measure can be found, the safety factor can be used as a substitute for the reliability index to monitor the safety of a specific dam.
In this paper, we propose measuring the ratio R = β/FS for different operational conditions over time, with other problem parameters remaining unchanged (such as nominal values and probability distributions of geotechnical parameters). After ratio R is computed at a minimum number of points in time, a confidence interval is evaluated for this ratio: where E() is the expected value operator. The lower and upper bounds are evaluated as: where Var() is the variance operator, and parameter k yields the desired confidence interval. For k = 2, for instance, one has a confidence of 95.5% that the actual ratio is contained within the bounds. Now assume that the mean ratio E(R) and the bounds in Eq. (8) have been evaluated, for some normal operational conditions, and the safety factor FS current is evaluated for another normal operational condition. The bounds on the estimated reliability index and the mean estimated reliability index are obtained as: Surely, the ratios and bounds in Eqs. (7) and (8) change for different dams, and for different probabilistic characterization of the seepage and stability properties of the dam. The relations above are assumed to be valid for a single dam in "stationary" normal operating conditions (NOC). Surely, the relations above are not valid in the presence of extreme loading events, such as rapid drawdown, extreme rainfall, degradation/failure of the drains or earthquakes.

Probabilistic modeling
The deterministic GeoStudio 2019 [19,20] software and the probabilistic StRAnD 1.07 [21] software were implemented using direct coupling to perform deterministic and probabilistic calculations. An earth dam in long-term steady-state and during a rapid drawdown was discussed by Siacara et al. [31,39]. In this research, an earth dam in normal operating conditions (NOC) was evaluated. The deterministic limit equilibrium method (LEM) (e.g. Morgenstern and Price) was used to find the factor of safety (FS). The reliability index (β) is evaluated in a specific time or in a time frame using structural reliability methods (e.g. FORM). The procedure implemented is described in this section. The flowchart shown by Siacara et al. [31,39], and the FORTRAN code were modified as shown in Fig. 1: Research Article SN Applied Sciences (2022) 4:99 | https://doi.org/10.1007/s42452-022-04980-7 (1) Initial data of the problem: studies of topography, geology, geotechnical (laboratory and field test), hydraulic and hydrological.  has a format ".gsz" which is a ZIP file that contains different input/output files. All input parameters are changed in the file extension ".xml", which is changed in every simulation or during a search for the design point (DP). The input parameters are organized/ changed in the original position of the ".xml" file for every simulation. For saving output results, file ".gsz" is compressed/uncompressed from the 7-zip software saving the FS to find the DP. All the output results (e.g., FS, PWP, phreatic surface, ..., etc.) are found in file extension ".csv" for every step of time. (11) Probabilistic results of the model: the STRAND_OUT-PUT.txt file of the StRAnD 2.00 software contains all the output results of the reliability analysis (number of evaluations and simulations, sensitivity coefficients at DP, reliability index, probability of failure, evaluation time, DP, ..., etc.). Visual results and additional information (e.g., pore water pressures in specific coordinates or critical slip surfaces) of the reliability results are found in the last ".gsz" file of the analysis. The probabilistic critical condition or minimum safety of the structure is represented by critical reliability index (β cr ) during the critical time (t cr-p ) (Fig. 2).
The FS is found from m to n steps of time or at a specific time of the transient analysis. From the FS, the StRAnD software performs the reliability analysis in both options of time (Fig. 2).

Earth dam of the study
The application of the probabilistic dam slope and flow analysis methodologies described earlier in the paper are illustrated in the analysis of a well-documented Brazilian dam. The dam is not identified due to privacy concerns. The reservoir filling was started in July 1970, and the construction of the dam was finished in August 1971. The dam has about 50 years of normal operation conditions (NOC), and it is classified as an old earth dam. The following study is an academic one; therefore, all the conclusions are only suggestions to the company in charge of the dam.
The dam was subject of extensive studies quantifying the soil properties from instrumentation, field and laboratory tests. However, we were not given access to measured data; hence, we could not evaluate the actual variability of the soil properties. These data were taken from the literature, but adjusted to the actual mean values provided in previous studies of the dam.
The earth dam is a dam of homogeneous type with different elements. The dam has a horizontal (approximately 1.8 m in height) and vertical filter (approximately 1.0 m in thickness) to reduce the pore water pressure within the dam, conducting the water downstream.
As the foundations presented good quality and favorable conditions for implementing the project, the concrete cut-off and the injection curtain initially foreseen in the project were not executed. Nevertheless, the dam counts on a small trench positioned parallel to its axis, excavated until the rock has some weathering to improve the waterproofing conditions of the structure. The dam foundation is of unweathered rock (granite and gneiss), and water level up to 4.0 m below the crest is the maximum water level. The saturated-unsaturated seepage and the stability of the dam slopes, considering the variability of soil parameters, is analyzed.
The plan view and cross-section are shown in Fig. 3 and critical cross-section A is presented in Fig. 4. Crosssection A is assumed to be the critical cross-section; it is used throughout the analysis herein. Different data and information on the earth dam studied are available [40,41]. The data are divided into three categories, according to type: field instrumentation, laboratory and field tests.
The case study advantage is that a reference system was installed during its construction, which allowed localizing the tests in space (X, Y and Z according to the three axes).

Dam condition of analysis
In this study, deterministic and probabilistic analyses are performed on a critical cross-section of the earth dam during NOC, where the embankment soil is saturated and unsaturated at different times of analysis. The deterministic parameters (nominal parameters) were quantified by field instrumentation, laboratory and field tests from previous studies of the dam, and the uncertainty in soil parameters was quantified from the literature [39].
The seepage analysis contains the representation of the SWRC and the hydraulic conductivity function predicted by using the Van Genuchten equation [29], and the Morgenstern and Price method [42,43] with the shear strength criteria of unsaturated soils [30] is used as LEM. The Van Genuchten fitting parameters (θ s ,θ r , a and n) were estimated using the Seep/W software database [44,45] for every dam material.  The dam body consists of silty-clays, the filter of sands, and the rip-rap and rock mass for protecting the fractured rock. For avoiding computational issues by differences in hydraulic conductivity, a transitional material of 5 cm was used between the dam body and filter. The random seepage analysis is characterized by the uncertainty of two SWRC fitting parameters (a, n), the saturated and residual volumetric water contents (θ s , θ r ), and the saturated hydraulic conductivity parameter (k s ). The random stability analysis is characterized by uncertain specific weight (γ), effective cohesion (c′), effective friction angle (ϕ′) and the angle that increases shear strength (ϕ b ). The uncertain parameters of the stability and seepage analyses are assumed to have normal (N) or lognormal (LN) distributions. The mean (µ) and the coefficient of variation (COV) of all parameters are shown in Table 2. In terms of COV of soil parameters, Phoon [16] reported different range of values and Siacara et al. [31,39] applied it in earth dams.
The SWRC and hydraulic conductivity curves corresponding to the mean values in Table 2, are shown in Fig. 5a, b, respectively.
The following remarks are made about the COV values and r k ratio in Table 2: 1. The COV of k s , θ s , θ r , a, n, γ, c′, ϕ′ and ϕ b for all the different soils variations was compilation of Siacara et al. [39]. 2. The r k ratio for existing dams varies from 1 to 15 [26,[46][47][48],Cruz [2] suggests a value of 1 for non-cohesive soils and from 4 to 10 for cohesive soils. In compacted earth fills, the r k ratio may exceed 20 [49]. Fell et al. [50] and UCACE [51] give r k ratios in the order of 1 to 100, covering the possible range of expected field conditions. A common value of 3 for r k ratio is suggested by Verbrugge and Schroeder [52]. For compacted soils, the r k ratio for non-cohesive soils varies from 1 to 40, and for cohesive soils from 0.4 to 4.1 [53]. From laboratory tests, the r k ratio varies from 1 to 4.1 [53]. Leroueil et al. [54] found r k ratios from 1 to 1.4 in compacted soils in laboratory, and Smith and Konrad [55] found an r k ratio of 5 using geostatistical analysis of construction control data from the core of an earth dam.
Concerning anisotropy, the data available did not allow using a random characterization of the anisotropy coefficients of the materials. The anisotropic coefficient (r k ) was considered a constant for the different materials. The r k ratio was assumed as 5 for dam body and transition, and 1 for filter and rock-mass, according to the literature on cohesive and non-cohesive soils ( Table 2).
The USBR [56] indicates that typical dams will have r k ratios ranging from 2 to 10, with higher values relating to higher water contents during placement. Older dams, such as those constructed in the early twentieth century or by hydraulic fill methods, may have anisotropy as high as 50 due to stratification during placement and earlier compaction methods that did not emphasize mixing and discing. However, coarse-grained materials, such as rockfill shells or filter and drain materials, are typically placed in thicker lifts without as much compactive effort, and they tend to have lower anisotropy. These types of soils are often assigned anisotropy values of 1.

Initial mean value analysis
Seepage and stability analyses were performed in the Geo-Studio software (Seep/W and Slope/W), for normal operating conditions (NOC), to obtain deterministic (mean value) results of pore water pressure (PWP) and factor of safety (FS) for the most critical slip surface. The mean value (µ) of seepage properties (k s , θ s , θ r , a and n) and properties involved in stability calculation (γ, c′, ϕ′ and ϕ b ), presented in Table 2, were used in this analysis. A two-dimensional analysis is performed herein. Table 3 shows the duration in days of the different analyses performed. The initial PWP condition is the key input information for transient seepage analysis in an earth dam, and it can be difficult to determine. The following methodology used herein will help to define the correct initial conditions. The first step (Table 3), an initial long term steadystate analysis was performed to determine the PWP for maximum water level conditions (845.0 m.a.s.l. or 75 m for numerical purpose), which is assumed at initial time (t = 0 days). The events described in Table 3 correspond to the actual normal operational conditions for the dam, following the official information available. The nomenclature of events states event duration in days (e.g., t = 570 days represents 570 days of analysis). This nomenclature corresponds to, but is independent of dam age. The height of the dam will be taken from the height of the numerical model in meters (m) and not from the meters above sea level (m.a.s.l.).
The initial seepage analysis (t = 0 days) yields the phreatic surface, PWP and total water head as shown in Fig. 6a. In this study, the dam stability is improved by negative  PWPs effects. The boundary between negative and positive PWPs (PWP = 0 kPa) is known as phreatic surface. The dam was modeled using triangular elements with of approximate absolute size of 5 m. The discretization of the mesh was according to the size of the element and the importance of the element in the analysis. A total of 33,825 nodes and 67,232 nodes were automatically created by the mesh generator available in Seep/W.
Taking into account flow measurements during construction, the infiltration trough the dam foundation was estimated at around 1.6 l/s, which, when increased by external contribution, has led to values of the order of 3.6 l/s. However, the real value was about three times higher than the one initially estimated. These differences eventually impacted the performance of the horizontal filter. The water level was found to be above the horizontal filter, and the downstream PWP was very different. This was verified with historical field measurements [57,58].
The numerical model employed herein was calibrated taking into account field measurements and instrumentation. The parameter used to calibrate the pore water pressures (PWP) is the saturated hydraulic conductivity (k s ) of the four materials of the dam. Real measurements of PWP are compared with the PWP of the numerical seepage analysis. The comparison was performed over time (150 continuous days) using different field instrumentations (different locations) inside the dam. The calibration was realized to find the most similar result between the field measurements and numerical results. The calibrated k s values were used in the numerical analysis, which yields representative results of real PWPs.
The second step (Table 3), the transient analysis (equilibrium time) is the variation of the reservoir (oscillation between 845.0 m.a.s.l. or 75 m and 827.0 m.a.s.l. or 57 m) over four years (from t = 0 days to t = 1460 days) to find the PWP equilibrium of the dam. The oscillation of the water reservoir is a variation from hydrological conditions, and by user demand request.
The behavior of the reservoir changes every day (there are always some water level variations), and the dissipation time of PWP is insufficient. In a transient analysis, the seepage analysis is performed to define the phreatic surface, PWP and total water head in every time step (in this case, one day). Figure 6b shows the seepage results for t = 500 days.
The third step (Table 3), the transient analysis is the variation in the reservoir (oscillation between 844.0 m.a.s.l. or 74 m and 831.0 m.a.s.l. or 61 m) over one year (from t = 1460 days to t = 1825 days). During this step, the earth dam is considered to be in equilibrium, and PWP values are more realistic. Figure 6c shows the seepage results for t = 1600 days. The stability of the dam is calculated from the downstream slope using the results of the three seepage steps mentioned above. The LEM is used for all stability calculations following the Morgenstern and Price's procedure. In every step of time, the seepage results (phreatic surface, PWP and total water head) are used to find the most critical slip surface. This is automatically found by the entry and exit specification technique (an indepth explanation is found in GeoStudio [19,20]. An initial configuration is defined by an 1.2 m discretization of the entry and exit slip surfaces, and radius tangent lines in the slip direction. The finite element mesh of seepage analysis, deterministic critical slip surface and the minimum factors of safety (FS) of the Slope/W software results are shown in Fig. 6a-c for times t = 0, t = 500 and t = 1600 days, respectively.
The variation of the water reservoir level as function of time produces different FS and critical surfaces for every time step, as shown in Fig. 7a, b. In Seep/W, an interpolated PWP was used to define the phreatic surface. All the slip surfaces found in Slope/W must be contained in the domain of the numerical model. If a slip surface follows the boundary, the strength parameters are taken from the materials overlying the base of individual slices. The initial condition of stability (t 0 = 0 days) yielded the initial safety factor as FS 0 = 2.13. During the equilibrium analysis (from t = 0 to t = 1460 days), the oscillation of FS is between FS max = 2.34 and FS cr = 2.13. During transient analysis time (from t = 1460 to t = 1825 days), the oscillation of FS is between FS max = 2.30 and FS cr = 2.18. Although the maximum water level is reached several times, the critical FS is found for t = 0 days, because the rise and fall velocity of the reservoir during NOC is not very large.
The mechanics behind the transient seepage analyses were developed more than 80 years ago. However, it is difficult to put the mechanics into immediate practical use owing to the large number of numerical calculations required. In recent years, transient seepage analysis programs, such as SEEP2D1, SEEP/W, and SLIDE, allowed these analyses to be conducted by personal computers in relatively modest execution times [59]. In this study, we consider unsaturated soils and employ long term steady-state analysis as an initial condition, followed by four years of equilibrium analysis, to find more realistic PWP before the actual transient analysis.
A minimum FS = 1.5 for long-term steady-state analysis and FS = 1.1 to 1.3 for rapid drawdown analysis in earth dams are recommend by different authors [2,50,60] and specialized organizations [61][62][63][64][65]. In this study, all the FSs of deterministic analyses meet these minimum stability criterions.

Reliability analyses and results
Reliability analyses for dam equilibrium were performed at the four time periods identified as A, B, C and D Fig. 7 and Table 3. (between t = 1460 and t = 1825 days). The direct coupling (DC) between GeoStudio/StRAnD softwares was used to perform reliability analysis. The FORM method was performed using the HLRF algorithm to search for the design point (DP).
In the reliability analysis, the random seepage and stability properties (see Table 2) were considered for every material of the dam. The results are presented as follows: (1) Preliminary screening of random variable importance in NOC; (2) Normal operating conditions at time periods A, B, C and D; (3) Differences of pore water pressures; (4) Differences of critical surfaces; (5) Sensitivity of the random variables w.r.t. seepage and equilibrium analyses.

Preliminary screening of random variable importance in NOC
An advantage of FORM is the possibility of carrying out a sensitivity analysis through the direction cosines (α 2 ) at the Design Point (DP). The relative contribution of each random variable used in a reliability analysis is measured by α 2 . Large and low α 2 values represent the most important and irrelevant variables, respectively. Initially, a reliability analysis was performed with 34 variables (Table 2) of the earth dam, for four different times. These results (not shown here) reveal that, out of 34 initial random variables, only four dam body parameters (k sat , γ, c′ and ϕ′) and two filter parameters (k sat and ϕ′) are important, with individual contributions to reliability analysis greater than 0.1%. The other 28 random variables resulted in nearly zero sensitivity coefficients (α 2 ≈0); hence, these random variables have negligible contributions to the computed failure probabilities and were considersed deterministic in the remaining analyses.
The numerical reliability analyses were performed by direct coupling, using a dual core workstation computer, with processor speed of 2.

Analysis for normal operating conditions (NOC)
During normal operating condition (NOC), we can obtain both deterministic FS and probabilistic β of the unsaturated earth dam. Figure 7 shows the variation in water levels over the analysis time, as well as the safety factors found in the deterministic (mean value) analysis. The reliability analysis is performed assuming the reservoir water level known at the end of each day; the end condition of one day is the starting condition for the next day. In dam practice, water levels at the end of any day are well known, by way of manual measurements (ruler inserted into water reservoir), hydraulic installations (capacity of spillways) or water demand (electricity generation or water supply). The reliability analysis is performed using the FORM method in the periods A (from t = 1466 to t = 1481 days), B (from t = 1600 to t = 1615 days), C (from t = 1645 to t = 1660 days) and D (from t = 1810 to t = 1825 days), as shown in Figs. 10 and 11. The four time periods of analysis (A to D) have lengths of 15 days (Table 3). In NOC analysis, the FS and β follow similar trends for the same step of time (Fig. 8). Time periods A, B and D show a negative tendency of the P f curve, but period C has a positive tendency. The FS and β curves have an opposite behavior to the P f curves for the four periods of analysis, as expected. The determination of FS, β and P f for every day of analysis reduces uncertainty regarding the behavior of the dam (e.g., t = 1470, 1610, 1655 and 1820 days), as shown in Fig. 8. Figure 8d shows two results for specific times (t = 1820 and 1825 days) with a difference of 5 days and a 0.66 m drop of the reservoir. The difference between results are these two points in time, ΔFS = 4.3 × 10 -3 , Δβ = 2.33 × 10 -2 and ΔP f = 4.17 × 10 -7 . These variations are important to prevent occasional unsafe behaviors.
In this manuscript, rapid drawdown is not considered part of normal operating conditions. Although rapid drawdown is typically considered the most critical equilibrium situation, the normal operating conditions could present other critical situations. An example is pore water pressure increases in the old dam by continuous rainfall and/or decrease of the efficiency of the dam filter. These situations lead to critical conditions of the dam, which can be detected using constant monitoring, or can be predicted using reliability analysis, as demonstrated herein.
The methodology developed herein targets the long term ongoing safety analysis of dams, which is based on combining transient seepage analysis with equilibrium analysis, and in relating safety factor changes with changes in failure probabilities. Figure 8 shows how FS, β and P f curves change during the normal operating conditions. As stated in Sect. 2.9, it may be possible to estimate reliability index β from the observed ratios R = β/FS of for this dam and for the random parameters considered in Table 2. Based on the results in Fig. 8, and using Eqs. (7) and (8), the ratios were found as R upper = 2.12, E(R) = 2.00 and R lower = 1.88. Based on these ratios, and on the FS values computed for the whole analysis period, the reliability index bounds were computed. The results are first compared for the four time periods (A, B, C and D) for which actual β values are known (Fig. 9). In Fig. 9, actual β values are observed to indeed be within the bounds given by Eqs. (7) and (8). This could be expected, since the ratios above were computed in the same time periods (A, B, C and D). The usefulness of the proposed bounds can be appreciated when observing the whole (deterministic) analysis time, in Fig. 10. In this figure, the bounds are employed to estimate the maximum and minimum values of β for the whole interval, including the times between time periods A, B, C and D. From these bounds, it is estimated that the maximum and minimum (critical) values of β are β max = 4.87 at t = 1626 days, and β cr = 4.11 at t = 1735 days. Hence, this initial screening can be used to perform a reliability analysis at the most critical section, around time t = 1735 days. This was done herein, after identifying t = 1735 as the critical time. The actual reliability analysis at time t = 1735 resulted in β real = 4.07; a difference of only Δβ = 0.04 to the estimated value. This confirms, by observation, that the bounds proposed in Eqs. (7) and (8) can be used for screening the most critical section in NOC. Note that t = 1735 is critical because it is the end of a longer period for which the reservoir had been full or nearly full; hence, pore water pressures were high in a significant part of the dam.
Considering the dam performance levels listed in Table 1 [36] and the reliability index values computed herein, the analyzed dam is classified as between Good and High. This includes the estimated critical value β cr . Equations (7) and (8) are alternatives to estimate β from the FS as is shown in Fig. 10. Although this approach helps geotechnical engineers to estimate reliability values from deterministic results, these equations not avoid to perform reliability analysis. The observed ratios R = β/FS of every dam can be improved in accuracy with more results from the reliability analysis, and a calibration with a high number of results gives more reliable estimation. These ratios R can be defined during the first years of NOC of the dam, and it can help during the safety monitoring of a dam. The principal advantage of using these ratios R are performing deterministic transient analysis to estimate the reliability results in function of time, avoiding more computational efforts. During NOC, the estimation of β is helpful (e.g. a dam during large increases or decreases of the reservoir; the state of the dam is known immediately). Factors of safety and reliability index complement each other and hence both are more useful than knowing either one alone.

Differences of pore water pressures
Differences between deterministic and probabilistic analyses in terms of seepage and stability are further investigated. In every step of time t k , different PWPs, phreatic and failure surfaces are illustrated over the dam. The results of the analyses (deterministic and probabilistic) are illustrated in the GeoStudio software. Deterministic results are represented by the mean values of the geotechnical variables (Table 2) and probabilistic results are represented by the DP values of the same geotechnical variables.
In Fig. 11, the deterministic and probabilistic PWPs of the seepage analysis are compared, for times t = 1604 and t = 1656 days. The critical equilibrium situation occurs for higher PWPs; therefore, the probabilistic PWPs in the DP are higher than the deterministic PWPs. The variation in PWPs depends on the time period of analysis, and on the velocity of water level change. The differences between results in Fig. 14a, b are mainly explained by use of mean value of k s (1.0 × 10 -6 m/s) and the design point value of k s (0.84 × 10 -6 m/s for t = 1604 and 1.24 × 10 -6 m/s for t = 1656 days). This difference is estimated to be larger in time period C.
Two cross-sections identified in Fig. 11 (A-A′ and B-B′) are presented in Fig. 12 to compare the PWPs at times t = 1604 and t = 1656 days. The points T and U identified in Fig. 11 are represented in Fig. 13, where the PWPs are shown in every time step of time periods B and C. Significant differences are observed between deterministic and probabilistic PWPs in time periods B and C over time. Hence, the differences between deterministic and probabilistic analyses can be explained by the effects of seepage in the days preceding the equilibrium analysis, and dissipation of PWPs in water level falls.

Differences of critical surfaces
A comparison of the corresponding critical deterministic and probabilistic slip surfaces, for times t = 1604 and t = 1656 days, is presented in Fig. 14a, b, respectively. The geometry of the deterministic critical slip surface does not change significantly as function of time, as seen in Fig. 6. The critical deterministic slip surface is located between the rock foundation and the rock mass downstream.
The reliability analysis using FORM involves a search for the design point (DP). When the random variables assume the values corresponding to the DP, the slip surface with the highest probability of occurrence is obtained at every time step. This is called the probabilistic slip surface.
As observed in Figs. 14 and 15, there are subtle differences between the PWPs and phreatic surfaces observed in the deterministic and probabilistic analyses. As discussed in the next session, this is a consequence of the (2022) 4:99 | https://doi.org/10.1007/s42452-022-04980-7 Research Article small contribution of random seepage parameters to the reliability problem. However, there are significant differences between deterministic and probabilistic slip surfaces, as observed in Fig. 14. The probabilistic slip surface is located above the downstream rock mass, in contrast to the deterministic slip surface. In transient NOC analyses, these results suggest that geometric differences between deterministic and probabilistic slip surfaces in homogeneous earth dams are much larger than differences observed in long-term steady-state analyses [31,39,66,67]. These differences need to be studied for zoned dams where the properties of the materials change in function of the geotechnical design. In normal operating conditions (NOC), safety factors and reliability indexes at a given time depend on the history of the reservoir water level and PWP conditions. Finding the real condition of earth dams is hard work in a deterministic approach, and using of probabilistic approaches is a challenge for geotechnical engineers. This helps to see the importance of performing reliability analyses during the design and for assessing the safety of old dams. Figure 15 illustrates sensitivity coefficients for the six most relevant random variables, under NOC and for time periods A, B, C and D. Firstly, the sensitivities are observed not to change significantly in time, which is quite different from the case or rapid drawdown [31]. This also adds to the argument that five days of previous seepage analysis  is sufficient for equilibrium analysis in NOC, whereas the whole time interval had to be considered in the rapid drawdown analysis [31].

Sensitivity of random variables
The results in Fig. 15 also show that, of the six random variables with the greatest contribution to failure probabilities, the dam body friction angle ϕ′ is the most important for all the NOC time periods considered herein. The second most important variable is dam body effective cohesion  c′. The other two parameters of the dam body (γ and k s ) and two parameters of the filter (k s and ϕ′) have smaller importance (α 2 ≠ 0). Studies on stochastic pore pressure variability [68], transient analysis of rapid drawdown [31] and stochastic hydraulic conductivity [69], as well as other studies [70,71], found seepage properties to have great importance to reliability analyses. The reliability index decreases as the COV of conductivity hydraulic increases. A comparison of our results with the references above reveals that the importance of the seepage parameters (K s , θ s , θ r , a and n) of the materials in our study are smaller or nearly zero (α 2 ≈0).
In this case, the geometry of the dam, the geometry and position of the filter, the phreatic surface and critical slip surface produce seepage properties which have smaller or nearly zero (α 2 ≈0) incidence. The critical slip is located downstream of the dam, and the importance of phreatic surface is less relevant when the reservoir is at lower water level.
In deterministic analyses, it is well known that the increase of phreatic surfaces downstream of the dam reduce safety factors. In the NOC reliability analysis, this results in reduction of β's and in a negative α 2 value for dam body k s . Water flow through the lower part of the critical deterministic and probabilistic surfaces (Fig. 14a, b) helps to explain the smaller relevance of random seepage parameters in the reliability analysis.

Concluding remarks
In this paper, reliability analysis of an existing old earth dam was performed for normal operating conditions (NOC), considering transient seepage conditions. A direct coupling (DC) between the deterministic Geo-Studio 2018 package (Seepage/W and Slope/W software) and the structural reliability StRAnD software was employed. Numerical analysis of seepage and stability were performed. The first-order reliability method (FORM) method was used to find the design point in reliability analyses.
For new dams, it is numerically possible to model all the history of the dam, in deterministic and probabilistic approaches, with admissible computational cost. However, an old dam is analyzed with a different approach. To find the most realistic pore water pressure (PWP), a long-term steady state analysis is performed to find the initial condition, followed by four years of equilibrium analysis before each transient seepage analysis. These PWPs are calibrated based on dam instrumentation readings. A computational analysis from the empty dam to the actual state of the dam is unfeasible, and reliability analysis would be very hard using FORM. In this study, five days of previous seepage analysis were found to be sufficient, starting from the mean value condition, for performing reliability analysis of a dam under NOC.
Identifying the critical time corresponding to minimum reliability is very costly. A simple empirical equation was proposed comparing safety factors (FS) and reliability indexes (β) for a single dam in NOC. This equation yields the expected (minimum and maximum) reliability indexes, based on FS calculations for the whole life of the dam. Estimated minimum FS values can be used to find the (approximated) critical time, at which actual reliability analyses should be performed. Note also that the critical time corresponding to minimum FS can be slightly different from the critical time corresponding to minimal β [31]. The empirical relation is valid for a single dam, and for a stationary description of random variables.
The initial analysis listed 34 random geotechnical parameters: five seepage and four equilibrium parameters for each material. Each day of random seepage analysis with 34 parameters took about 20 h to compute. Sensitivity analyses revealed the six parameters with the greatest impact on evaluated failure probabilities; a reduction in random variable dimensionality reduced processing time for a daily transient analysis to about 6 h. Using deterministic (mean values) data of laboratory and field test, and statistical data from the literature, the transient analysis reveals that four dam body parameters (k sat , γ, c′ and ϕ′) and two filter parameters (k sat and ϕ′) presented the largest contribution to reliability analysis. The friction angle (ϕ′) has the greatest impact on the reliability analysis for equilibrium. Different conditions of the dam (e.g., the geometry of the dam and position of the filter, the phreatic surface and critical slip surface) produce random seepage properties with smaller or nearly zero (α2 ≈ 0) contribution. The critical slip is located downstream of the dam, and the importance of phreatic surface is less relevant when the reservoir is at lower water level.
The cumulative effect of random k sat (mainly) in transient analysis produces worse critical seepage results (phreatic surfaces, PWP and total water head) and stability results (critical slip surface) in probabilistic analysis than in deterministic (mean value) analysis.
The most probable slip surface, found in the reliability analysis, is not the same as the critical slip surface found in a mean value analysis for the minimum factor of safety (FS). The difference is due to seepage parameters and the geometry of the dam.
Considering target reliability indexes suggested in ICOLD [72], the expected performance level of the studied dam is Good.
A limitation of the study was the determination of different conclusions using only four time periods. The Data availability Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.