The Reliability of Stored Water Behind Dams Using the Multi-Component Stress-Strength System

: Dams are essential infrastructure for managing water resources and providing entry to clean water for human needs. However, the construction and maintenance of dams require careful consideration of their reliability and safety, speciﬁcally in the event of extreme weather conditions such as heavy rainfall or ﬂooding. In this study, the stress-strength model provides a useful framework for evaluating the reliability of dams and their ability to cope with external stresses such as water pressure, earthquake activity, and erosion. The Shasta reservoir in the United States is a prime example of a dam that requires regular assessment of its reliability to guarantee the safety of communities and infrastructure. The Gumbel Type II distribution has been suggested as a suitable model for ﬁtting the collected data on the stress and strength of the reservoir behind the Shasta dam. Both classical and Bayesian approaches have been used to estimate the reliability function under the multi-component stress-strength model, and Monte Carlo simulation has been employed for parameter estimation. In addition, some measures of goodness-of-ﬁt are employed to examine the suitability of the suggested model.


Introduction
In reliability analysis, the study of the failure time of a component or a system is of great concern.Decision-makers in industry or governments encourage continuity in the enhancement of systems reliability since it measures the functioning of systems and predicts their outcomes in a better way.One of the main problems is estimating the parameters of the stress-strength function R = P(Y < X), where Y is the stress and X is the strength, and both are considered random variables.Components or systems are subject to failure if their strength is less than the stress applied to them at any time.
A multi-component system, with k-independent and identically distributed (iid) strength components and one common stress, is in a functional state if at least s components are surviving at the same time, where 1 ≤ s ≤ k, and this is known as an s-out-of-k: G system.Two examples of multi-component systems are: first, the suspension bridge is suspended by a k-vertical cable pair, this bridge will survive if at least s-vertical cable pairs are functioning, i.e., not damaged.Second, the engine with six cylinders will work properly if at least four out of six cylinders are functioning.
Estimating the reliability of a system with the multi-component stress-strength model is of great importance and has many applications, see for example Ahmad et al. [22], who used the Power Lomax distribution under a progressive censoring scheme.Johari et al. [23] studied the reliability analysis of ground soil layers using the cross-correlation method.More reliability inference studies are presented in [24,25].
In this work, we consider the stress-strength model applied to reservoirs or artificial lakes behind dams.Dams are very important for providing water for different purposes during droughts such as drinking, cooking, washing, and irrigation, as well as producing hydroelectric power.They help in storing water in reservoirs during excess water flow and release it during times of low flow.Besides irrigation and other consumption purposes, water stored in these reservoirs is subject to other loss factors such as evaporation.We consider a model that is related to excessive drought times, and it is also claimed that in a specific region, if the stored water in a reservoir in August of the previous year is less than the water amount in the reservoir for at least two of the next four years, then no excessive drought will occur.In this model, X represents the strength and Y represents the stress.
In this study, we aim to estimate the stress-strength function's parameters R s,k which follow the Gumbel Type II model.This model was invented in 1958 by the German mathematician Emil Gumbel  and was useful for predicting the likelihood of climatic events such as annual flood flows (which is the case of this study), earthquakes, and other natural disasters.According to Abbas et al. [26], it has also been demonstrated to be sufficient for describing the component's life expectancy.In their book, Kotz and Nadarajah [27] investigated the Gumbel distribution with the aim of applying it to the analysis of datasets ranging from wind speed to flood data.
Several authors focused on the Gumbel Type II distribution in their research for its wide range of applications.For example, Nadarajah and Kotz [28] studied the Gumbel distribution's main properties.Feroze and Muhammad [29] performed Bayesian inference for the Gumbel Type II distribution under twice censored samples with different loss functions.Mansour and Aboshady [30] discussed different estimation methods using the Gumbel Type II distribution under the hybrid Type II censored scheme, they applied this model to study the performance of the insulating fluids for the breakdown voltages.
One of the basic motivations for using the Gumbel Type II distribution for modelling real data is its ability to keep the properties of a wide range of distributions, including the Weibull, Fréchet, and Gumbel distributions.This makes the Gumbel Type II distribution a flexible alternative for modelling data that are not compatible with some distributional forms.Another motivation for using the Gumbel Type II distribution is its ability to model both heavy-tailed and light-tailed distributions, as well as distributions that obtain positive or negative skewness.Moreover, the Gumbel Type II distribution is often used to model extreme events such as floods, hurricanes, and earthquakes, where accurate modelling of the distribution is essential for risk assessment and management.The Gumbel Type II distribution has been found to support a good fit for many extreme events, making it a preferable choice for modelling such events.
Even though several types of research have been performed on the statistical inference of the Gumbel Type II distribution, there is still a persistent need for more work in the field of multi-component stress-strength systems with Gumbel Type II distribution.In this study, we work on two goals: first, to find classical and Bayesian estimators of multi-component stress-strength functions and assess the estimators by providing a simulation study and applying suitable numerical techniques.Second, apply the new model to fit real data collected from reservoirs behind the Shasta dam.
The procedure for organizing the process of this work is schematically shown in the flowchart in Figure 1.

Reliability Function
In this section, a closed form is obtained for the reliability function for a stress-strength model.Assuming the probability density function (PDF) and the CDF for the Gumbel Type II distribution, with a random variable X defined as, and and given the reliability function for a stress-strength model, that was defined by Bhattacharyya and Johnson [31], as where X 1 , X 2 , . . ., X k are iid random variables that follow the Gumbel Type II distribution with parameters α 1 and β 1 with cumulative distribution function (CDF) F(X), and are affected by random stresses Y that follow the Gumbel Type II distribution with parameters α 2 and β 2 having CDF G(Y).Substituting by the CDF in Equation (2) into R s,k defined in Equation (3), we obtain Let u = x −α 2 and substitute it into Equation (4) to obtain Equation ( 5) By using the binomial expansion for [1 Since the Maclaurin expansion for [e −β 1 u , then Equation ( 6) is rewritten as Assuming z = β 2 u, Equation ( 7) is reduced to

Maximum-Likelihood Estimation
This section is devoted to deriving the maximum likelihood estimators (MLEs) for the reliability function's R s,k parameters.The basis for finding the estimators of the parameters is the log-likelihood functions, given the data.The MLEs have been used by many authors to derive the estimators of the parameters due to several advantages such as simplicity, unbiased for larger samples, and acquiring smaller variance.Furthermore, it can be developed for a large variety of other estimation methods.For more information on the likelihood theory, see Azzalini [32].
To find the MLEs of R s,k , we start by obtaining the MLEs for the parameters α 1 , α 2 , β 1 and β 2 .In this model, samples can be constructed as

Observed strength variables Observed stress variables
Hence, the likelihood function for the observations can be written as where Equation (9) can then be written as From Equation (10), the log-likelihood function can be derived as follows Computing the first partial derivatives for with respect to α 1 , α 2 , β 1 and β 2 and equating with zero will give the four equations below nk nk n Straightforward from Equations ( 14) and ( 15), the estimators for β 1 and β 2 can be given as and β2 = n Dragging the two estimates in Equation ( 16) into Equations ( 12) and ( 13), respectively, and solving will give estimates for α 1 and α 2 , α1 and α2 , respectively, and substituting by α1 , α2 , β1 and β2 in Equation ( 8), we obtain the MLE Rs,k as follows

Fisher Information Matrix
In this sub-section, the asymptotic confidence interval (ACI) for the reliability function will be derived using the Fisher information (FI) matrix.The concept of the FI matrix is based on the missing value principle which was introduced by Louis [33] and is defined as follows: Observed information = Complete information − Missing information.
The asymptotic variance-covariance of the MLEs α1 , α2 , β1 and β2 are derived from the entries of the inverse matrix of the FI matrix Unfortunately, obtaining an exact closed form for the previous expectations is very complicated.Hence, the observed FI matrix Îij = − ∂ 2 (Φ)/∂φ i ∂φ j Φ= Φ.This is obtained by dropping the expectation op- erator E and using it to construct the confidence intervals for the unknown parameters.
The entries of the observed FI matrix are the second partial derivatives of the log-likelihood function, which is easily obtained.Therefore, the observed FI matrix is given by , (18) By inverting the information matrix Î(α 1 , α 2 , β 1 , β 2 ), the approximate asymptotic variance-covariance matrix [ V] for the MLEs can be obtained as: Assuming some regularity conditions, (α 1 , α2 , β1 , β2 ) will be approximately distributed as a multivariate normal distribution with mean (α 1 , α 2 , β 1 , β 2 ) and covariance matrix [34].Then, the 100(1 − γ)% two-sided confidence intervals of α 1 , α 2 , β 1 and β 2 can be given by where Z γ 2 is the percentile of the standard normal distribution with a right-tail probability γ 2 .
To construct the ACIs of the reliability function, R s,k , it is necessary to compute its variance.The MLE of the R s,k is asymptotically normal with mean R s,k and its corresponding asymptotic variance is given as Then, the 100(1 − ϑ)% two-sided confidence intervals of R s,k can be given by

Bayesian Estimation
Another method for obtaining the estimates for the distribution parameters and the reliability function is discussed in this section and is known as the Bayesian estimation.Before collecting and organizing the data, the joint prior distribution should be assumed, and what distinguishes this method is that the prior knowledge is merged in the solution steps.The Bayesian estimates for the four parameters α 1 , α 2 , β 1 and β 2 in addition to the reliability function R s,k is obtained under the squared error loss (SEL) function.First, the prior knowledge of the parameters α 1 , α 2 , β 1 and β 2 is assumed to follow gamma distribution as follows assuming that all the hyper-parameters a i and b i , i = 1, 2, 3, 4 are known and non-negative.
It can be noted that one reason for choosing this prior density is that Gamma prior is flexible in its nature with a non-informative domain, especially if the values of the hyperparameters are assumed to be zero, for more details on selecting priors one may refer to Kundu and Howlader [35], Dey and Dey [36], and Dey [14].
Using the likelihood function in Equation ( 10) and the prior distribution for the parameters α 1 , α 2 , β 1 and β 2 assumed in the previous equations, the posterior distribution, denoted by π * (α 1 , α 2 , β 1 , β 2 | x ¯, y ¯), for these parameters can be derived as follows Given a parameter φ which is estimated by φ, the symmetric loss function SEL function assigns equal losses for both over-and under-estimations, which can be defined as As a result, the Bayes estimate g(α 1 , α 2 , β 1 , β 2 ) under the SEL function can be writ- ten as ĝBS α 1 , α 2 , β 1 , , where The joint posterior density function of α 1 , α 2 , β 1 and β 2 can be obtained as follows The Bayesian estimate of R s,k , under the SEL function, is the mean of the posterior function in Equation (25) and can be written as It is clear that the integral in Equation ( 27) is difficult to be calculated analytically.Therefore, the Gibbs sampling method is used to obtain the Bayesian estimator for the reliability function R s,k .

Gibbs Sampling
The Gibbs sampling method is a special case of the Monte Carlo Markov Chain (MCMC) and can be used to perform the Bayes estimate of R s,k numerically, in addition to the its related credible interval (CRI).The key idea in Gibbs sampling is to generate samples for the required parameters from the posterior conditional density function, given in Equation (26).Then, the posterior conditional density functions of α 1 , α 2 , β 1 and β 2 are given as and respectively.It is difficult to obtain the conditional density function of α 1 , α 2 , β 1 and β 2 .Therefore, the Metropolis-Hasting (M-H) algorithm, proposed by Metropolis et al.
[37], is applied using the normal proposal distribution for generating random samples from the posterior density of α 1 , α 2 , β 1 and β 2 .The steps of Gibbs sampling are described as follows: 1.
Start with initial guess α

3.
Using the following M-H algorithm, generate α 1 and β (l) with the normal proposal distributions where Var(α 1 ), Var(α 2 ), Var(β 1 ) and Var(β 2 ) can be obtained from the main diago- nal in the inverse Fisher information matrix.
To compute the CRIs of α 1 , α 2 , β 1 , β 2 and R s,k , ψ The first M simulated variants are discarded in order to ensure convergence and remove the affection of initial value selection.Then the selected samples are ψ (i) k , j = M + 1, . . ., N, for sufficiently large N.
Based on the SEL function, the approximate Bayes estimates of ψ k is given by

Real Data Analysis
In this section, the reliability function is estimated using the MLE and Bayesian estimation methods, where the data under consideration are obtained for the water capacity in the Shasta reservoir in the United States.The view of Shasta Lake during the season of floods in addition to a general view of the Shasta dam are shown in Figure 2. To consider the scenario of the excessive drought, we will focus on the total amount of water in the period from 1980 until 2019.Our claim is that an excessive drought occurs if the total amount of water in August in two years of the next four years is less than the amount of water filling the reservoir in December of the preceding years, otherwise, no excessive drought will occur.This problem was previously studied in different contexts, i.e., see Fatma [19] and Akram [16].The source of the data is available in [38].
For computational simplicity, the water amount in the reservoir for any given month is divided by the total capacity of the reservoir and the data will then be as follows: The first row of the observations in matrix X represents the storage amount of water (divided by the total capacity storage of the lake) in August 1980-1983, respectively, while the second row represents this in August 1985August -1988, respectively, and so on, with a sample size of 32 observations.Whereas the amount of water in December 1984, 1989, up to 2019 is represented by matrix Y, with a size of 8 observations.It was found that both X and Y values are well-fitted to the Gumbel Type II distribution.The Gumbel Type II distribution was chosen primarily because of its suitability for predicting the possibility of climatic events such as annual flood flows (which is the case of this study), earthquakes, and other natural disasters.However, fitting outcomes for certain classical models, such as the log-normal distribution, are poor, with p-values of 1.57936 × 10 −11 and 0.00194213 for the X and Y datasets, respectively.While for the Gumbel Type II distribution, the Kolmogorov-Smirnov (KS) distance for X with the estimated parameters is 0.178521 and the corresponding pvalue is 0.230418.Furthermore, for Y, the KS distance with the estimated parameters is 0.263293 and the corresponding p-value is 0.550885.A comparison between the empirical distribution of the dataset and the survival function of the Gumbel Type II distribution is presented in Figure 3.For the complete dataset, the MLEs for the parameters, α 1 , α 2 , β 1 , and β 2 , and the reliability function, in addition to the Bayesian estimation with respect to SEL function are displayed in Table 1.From Table 1, the estimated values are relatively close to each other, which indicates the good performance of the estimators.In addition, the two estimators for the reliability function R 2,4 seem to be very close and approximately equal to zero, indicating no excessive droughts in these periods.The convergence for the estimated parameters using the MCMC method with 1000 iterations is shown in Figure 4.

Simulation Study
A Monte Carlo simulation analysis is carried out in this section comparing the performance of the MLE and Bayesian estimates of the distribution's parameters α 1 , β 1 , α 2 , and β 2 under various scenarios.The comparisons are made using the mean-squares-error (MSE) criterion.All computations are carried out in Mathematica 12, and all the results are obtained using 1000 Monte Carlo samples.We use various sample sizes (n) and (s, k) values in the simulation setup.Various parameter values are selected with different sample sizes (n) such as n = 30, 40, 50, 100, and 150 and α 1 , β 1 , α 2 , and β 2 are assumed to be (0.2, 0.5, 0.05, and 200) and (0.3, 0.4, 0.05, and 200), respectively.
Tables 2 and 3 show that in most cases, the Bayesian estimators of the parameters have a large bias compared with the MLEs for small sample sizes.However, when the sample sizes grow larger, all of the estimators exhibit small biases.In terms of the MSE criterion, it is clear that as the sample size increases, the MSEs for the estimates of α 1 , β 1 , and α 2 decrease, as expected.However, it was noticed that the MSE for β 2 does not decrease as the sample size increases; hence, sometimes some values of MSE violate this pattern due to some causes such as the numerical solution of a certain number of non-linear simultaneous equations.It should be noted that informative priors improve the performance of Bayesian estimates in a reasonable way.As sample sizes grow larger, the MSE values for all estimators become nearly identical.To emphasize the performance of the proposed methods, another criterion is applied which is the coverage probability.The coverage probability indicates how many times a confidence interval contains the initial value of the estimated parameters through the number of simulated samples.In addition, the average interval lengths for both the ACIs and CRIs are calculated and tabulated in Tables 4 and 5.It is noted that for various values of α 1 , β 1 , α 2 , and β 2 , for n, s, and k, the CRIs perform better than the ACIs through the coverage probability and the average lengths of intervals.

Conclusions
In this study, the stress-strength model describing the amount of water held in a certain reservoir over time is modeled parametrically using the Gumbel Type II distribution.It introduces a mathematical procedure for obtaining a closed form for R. In addition to the Bayesian technique, the maximum likelihood estimation for the parameter R is carried out.

Figure 1 .
Figure 1.A flowchart indicating the research process.

Figure 2 .
Figure 2. View of Shasta Lake during the season of floods and a plan view for the dam.

Figure 3 .
Figure 3. Empirical and fitted survival functions for the two datasets X and Y.