BMA is a statistical analysis method and can be used to infer a probabilistic prediction. It is a statistical analytical method that considers the uncertainty of the model itself. In this paper, the BMA method is applied to simulate the streamflow. Suppose that y is the forecasted variable, F = ( f 1 , f 2 , ⋯ f k ) is the all considered model predictions, and M = ( y 1 , y 2 , ⋯ y t ) is a given empirical probability distribution, where y t is the experience distribution probability at t time. Considering the inherent uncertainty of a model, the posterior distribution of y under a given empirical probability distribution condition can be presented as: p ( y | F , M ) = ∑ k = 1 k p ( f k | M ) p k ( y | f k , M ) (1) where p ( f k | M ) represents the posterior probability of a single distribution model f k , p k ( y | f k , M ) means the posterior distribution of y under a given model prediction f k and data set M.

Let ω k represent p ( f k | M ) with ∑ k = 1 K ω k = 1 . The posterior mean value and variance of the BMA prediction can be expressed as (Raftery et al., 2003; Duan et al., 2006): E ( y | M ) = ∑ k = 1 K p ( f k | M ) E [ p k ( y | f k , M ) ] = ∑ k = 1 K ω k f k (2) V a r ( y | M ) = ∑ k = 1 K ω k ( f k − ∑ i = 1 K ω i f i ) 2 + ∑ k = 1 K ω k σ k 2 (3) where σ k 2 represents the variable variance when the empirical frequency and model prediction are M and f k , respectively. The BMA mean value prediction can be calculated by the weighted average probability distributions for different optimal simulation. The two parts denoted by ∑ k = 1 K ω k ( f k − ∑ k = 1 K ω i f i ) 2 and ∑ k = 1 K ω k σ k 2 in Eq. 3 mean the error between models and error of the model itself individually. Compared to deterministic multi-model combinations, BMA possesses more reliability in reflecting the uncertainty of distribution models.

To effectively calculate the weight and variance of BMA, the expectation-maximization (EM) method is introduced in this paper. EM is an iterative calculation method and has been widely used (Bilmes, 1998) to calculate the maximum likelihood estimation, and the obtained results show good effect especially in dealing with a great number of missing data. In detail, the EM algorithm can be illustrated as follows.

Considering the stability and convenience of calculation, the EM algorithm uses log-likelihood function (Duan et al., 2006). Assume the prediction error in time and space is independent, the log-likelihood function can be formulated as follows (Duan et al., 2006): ℓ ( ω 1 , ⋯ , ω k , σ 2 ) = ∑ t = 1 T log ( ∑ k = 1 K ω k g k ( y t | f k , t ) ) (4)

Assume Z k , t is an unobserved variable, if at time t, f k is the best prediction model, then Z k , t = 1 ; otherwise, Z k , t = 0 . At any time t, only one of { Z 1 , t , ⋯ , Z k , t } is equal to 1 and the others are equal to zero. The key of applying the EM method is to switch steps between the expectation and maximization. It starts under an initial guess, such as assuming a value of θ ( 0 ) for parameter θ . In the expectation (E) step, Z k , t is estimated by the given current guess value of θ . In the maximization (M) step, θ is estimated by the given current values of Z k , t . These two steps are repeated until a predetermined accuracy meets the requirement. The detailed EM computation process is as follows.

Step 1. Initialization: S e t   j = 0 , ω k ( j ) = 1 k , σ k 2 ( j ) = 1 k ∑ t = 1 T ( ∑ k = 1 K ( y t − f k , t ) 2 ) T (5) where T is the number of data, and j is the number of iterations.

Step 2. Calculate the initial likelihood: ℓ ( θ ( j ) ) = log ( ∑ k = 1 K ω k ⋅ p k ( y | f k , M ) ) = log ( ∑ k = 1 K ω k ⋅ ∑ t = 1 T g ( y t | f k , t , σ k ( j ) ) ) (6)

Step 3. Implementation of the E-step operation: S e t j = j + 1 , Z k , t ( j ) = g ( y t | f k , t , σ k ( j − 1 ) ) ∑ k = 1 K g ( y t | f k , t , σ k ( j − 1 ) )   k = 1,2 , ⋯ , K , t = 1,2 , ⋯ , T (7)

Step 4. Perform the M-step operation:

Calculate the weights of the single prediction model: ω k ( j ) = 1 T ∑ t = 1 T Z k , t ( j ) (8)

Update the single forecast model variance: σ 2 ( j ) = ∑ t = 1 T Z k , t ( j ) ⋅ ( y t − f k , t ) 2 ∑ t = 1 T Z k , t ( j ) (9)

Update the value of the likelihood function ℓ ( θ ( j ) ) .

Step 5. Check convergence:

If ℓ ( θ ( j ) ) − ℓ ( θ ( j − 1 ) ) is less than or equal to 10 − 4 , then stop the operation; else return to the third step.

Readers can refer to McLachlan and Krishnan (McLachlan and Krishnan, 1997) for more detailed information of the EM algorithm.

In statistics, the kernel density estimation (KDE) is a non-parametric way to estimate the probability density function of a random variable. KDE is a fundamental data smoothing problem where inferences about the population are made, based on a finite data sample (Li et al., 2022). Among the many non-parametric methods currently used, the KDE method proposed by Guo et al. (1996) is the most widely used, and the effect is the most ideal.

The estimation formula of the univariate kernel probability density function: f ^ x ( x ) = 1 n h ∑ i = 1 n k ( x − x i h ) (10) where n is the length of observed data x i ; K ( ⋅ ) is the kernel density function; and h is the window width, it determines the variance of the kernel function.

Let u_i be the variable margin, θ be the parameter of the Copula function, the cumulative probability distribution of the three Archimedean Copula functions mentioned above can be expressed as follows:

(1) Gumbel-Hougaard Copula function

Correlation measure of random variables is used to describe the mutual dependence between random variables. There are many test metrics, in this study Kendall’s rank correlation coefficient τ and Spearman’s rank correlation coefficient ρ are used to evaluate the correlation of the streamflow variables. In detail, the corresponding Kendall’s rank correlation coefficient τ and Spearman’s rank correlation coefficient ρ s for Copula function C ( u , v ) can be expressed as follows: τ = 4 ∫ 0 1 ∫ 0 1 C ( u , v ) d C ( u , v ) − 1 (17) ρ s = 12 ∫ 0 1 ∫ 0 1 u v d C ( u , v ) − 3 = 12 ∫ 0 1 ∫ 0 1 C ( u , v ) d u d v − 3 (18)

To estimate the parameter of Copula function, the exact maximum likelihood (EML) method is used (Dupuis, 2007). If the joint distribution function of t-dimensional continuous random variables x 1 , x 2 , ⋯ x t can be expressed as C ( u 1 , u 2 , ⋯ u t ) , its maximum likelihood estimate can be calculated as follows (Dupuis, 2007):

Step 1. Establish the relevant likelihood function

The joint density function is expressed as: c [ ( u 1 , u 2 , ⋯ u t ) ; ( θ 1 , θ 2 , ⋯ , θ t ) , α ] = ∂ t C [ ( u 1 , u 2 , ⋯ u t ) ; ( θ 1 , θ 2 , ⋯ , θ t ) , α ] ∂ u 1 ∂ u 2 ⋯ ∂ u t (19) Likelihood function is expressed as: L [ ( θ 1 , θ 2 , ⋯ , θ t ) , α ] = ∏ i = 1 n c [ ( u i 1 , u i 2 , ⋯ , u i t ) ; ( θ 1 , θ 2 , ⋯ , θ t ) , α ] (20) Correspondingly, the log-likelihood function can be expressed as: ln ⁡ L [ ( θ 1 , θ 2 , ⋯ , θ t ) , α ] = ∑ i = 1 n ln ⁡ c [ ( u i 1 , u i 2 , ⋯ , u i t ) ; ( θ 1 , θ 2 , ⋯ , θ t ) , α ] (21)

Step 2. Solve the likelihood function ( θ ^ 1 , θ ^ 2 , ⋯ , θ ^ t ) , α ^ = arg max ⁡ ln ⁡ L [ ( θ 1 , θ 2 , ⋯ , θ t ) , α ] (22) where: ( θ 1 , θ 2 , ⋯ θ t ) are the unknown parameters of the distribution function for each margin; α is the associated unknown parameter in the Copula function.

In the procedure of hydrological frequency analysis, the monthly rainfall and runoff probability distributions in the Xiangxi River basin are first estimated by the Gamma, the generalized extreme value, and the lognormal distributions, respectively. And then, the BMA-based marginal distributions are obtained according to the three estimated distributions. Table 1 shows the fitting parameters of probability distributions.

Based on the weights and distribution parameters presented above, the BMA-based marginal distributions of monthly rainfall and runoff can be obtained. The Gamma distribution may account for a major proportion (99.99%) to produce the BMA-based marginal distribution of monthly rainfall; while the generalized extreme value distribution and lognormal distribution may account for almost the same proportion to produce the BMA-based marginal distribution of monthly runoff.

The comparison of empirical and generated marginal cumulative distribution functions (CDFs) for monthly rainfall and runoff is shown in Figure 6. It indicates that the BMA-based marginal distribution may appropriately represent the univariate rainfall and runoff probability distributions. In order to clearly clarify, the D, RMSE, and AIC values for the marginal distributions obtained by the four methods are also calculated and presented in Table 2. In addition, to compare with the non-parametric methods, the KDE method is also calculated in the same way. Results D show that only the GEV and BMA methods pass the K-S test (the upper boundary of D is 0.092 while alpha is 0.05) in both rainfall and runoff. The obtained results indicate the corresponding errors of the BMA method are relatively smaller suggesting the accuracy of the marginal distributions generated by the BMA method is very excellent.

After determining marginal distributions, the joint probability distributions of monthly rainfall and runoff in the Xiangxi River can be estimated by a Copula function. The estimation parameters for each Copula function are calculated based on the maximum likelihood estimation theory. In addition, according to the obtained parameters, the correlation coefficients can be calculated. Results are given in Table 3. It can be seen that the Kendall’s rank correlation coefficient ranges from 0.42 to 0.59 and the Spearman’s rank correlation coefficient ranges from 0.56 to 0.78. Therefore, it can be concluded that the monthly rainfall and runoff of the Xiangxi River have a relatively strong positive correlation.

According to Table 3, Gumbel-Hougaard-Copula-based joint distribution and Frank-Copula-based joint distribution are superior to the Clayton one based on the Kendall correlation and the Spearman correlation. Moreover, it shows that the estimation of the upper tail correlation coefficient should select the Gumbel Copula, the number is 0.6696. The estimation of the lower tail correlation coefficient should select the Gumbel Copula, the number is 0.7516. Obviously, there are both upper tail correlation and lower tail correlation in the Xiangxi River. This conforms to research by Yang et al. (2016). Goodness-of-fit tests of empirical joint CDFs and theoretical joint CDFs are calculated for further analysis (as shown in Figure 7).

Comparing the results shown in Figure 7, it can be known that the joint distribution by Frank-Copula is close to the one by Gumbel-Hougaard-Copula in quantifying the relevant characteristics of monthly rainfall and runoff of the Xiangxi River. The RMSE of the Gumbel-Hougaard-Copula is less than the Clayton copula. The Frank-Copula function has no tail correlation and cannot capture the tail correlation between variables according to Xue (2018). Gumbel Copula performs better in the upper correlation while Clayton Copula performs better in the lower correlation (Jondeau, 2016). Therefore, the Gumbel Copula function would be chosen to construct the joint distribution of monthly rainfall-runoff pairs. The corresponding results are plotted and shown in Figure 8. Figure 9 shows the combined probability density and return period of rainfall and runoff using the Gumbel Copula. It indicates that the extreme value of annual precipitation or annual runoff, the joint probability density is relatively small. When the annual runoff is constant, the greater the annual precipitation, the longer the time of return period is. In addition, the largest rainfall of observed data is 1,341.7 mm, the runoff of that year is 556.36 m³/s. The return period is about 10 years, which is reasonable.

In this study, the AD test and the CM test also have been selected to investigate the suitability of the BMAC-based joint distributions in describing the dependencies for different rainfall-runoff pairs. The results are displayed in Table 4, and statistics A n 2 are less than threshold values A n , 0.1 2 ( ∗ ) ( α = 0.1 ) and A n 2 ( ∗ ) ( α = 0.05 ) , where α is the significant level. Thus, the null hypothesis H 0 would be accepted. In total, it can be concluded that BMAC has a distinct superiority in modelling variable pairs.

In order to further clarify the efficiency of the BMAC method, a Maximum Entropy-Gumbel-Hougaard Copula (MEGHC) method proposed by Kong et al. (2015) has been applied for comparison. Accordingly, two joint distributions can be generated by the two methods, and the bias value between the two methods also can be obtained (shown in Figure 10). From Figure 10, the results are quite close to each other, and the great deviation value (i.e., absolute error) is 0.03, which illustrates that both of the two methods can be used to generate the joint distribution of the rainfall and runoff because of the best fitting effect. It also can be found that, under the condition of the CDF interval being 0–0.4, results obtained by the BMAC method are superior to the MEGHC method; while the CDF interval being 0.9–1, the MEGHC method may converge much faster. To some extent, the MEGHC method can better capture the characteristics of the upper tail dependence, which plays a great role in flood and drainage control and watershed design work (Kong et al., 2015). However, in the situation of representing the uncertainty of the model structure, it is hard to find a suitable distribution to capture the characteristic of the hydrologic variable. Comparatively, the BMAC method can obtain a synthetic simulation result especially in exactly reflecting the variables’ correlations.