Frequentist and Bayesian Approaches in Modeling and Prediction of Extreme Rainfall Series: A Case Study from Southern Highlands Region of Tanzania

Tis study focuses on modeling and predicting extreme rainfall based on data from the Southern Highlands region, the critical for rain-fed agriculture in Tanzania. Analyzing 31years of annual maximum rainfall data spanning from 1990 to 2020, the Generalized Extreme Value (GEV) model proved to be the best for modeling extreme rainfall in all stations. Tree estimation methods– L -moments, maximum likelihood estimation (MLE), and Bayesian Markov chain Monte Carlo (MCMC)–were employed to estimate GEV parameters and future return levels. Te Bayesian MCMC approach demonstrated superior performance by incorporating noninformative priors to ensure that the prior information had minimal infuence on the analysis, allowing the observed data to play a dominant role in shaping the posterior distribution. Furthermore, return levels for various future periods were estimated, providing guidance for food protection measures and infrastructure design. Trend analysis using p value, Kendall’s tau, and Sen’s slope indicated no statistically signifcant trends in rainfall patterns, although a weak positive trend in extreme rainfall events was observed, suggesting a gradual and modest increase over time. Overall, the study contributes valuable insights into extreme rainfall patterns and underscores the importance of L -moments in identifying the best ft distribution and Bayesian MCMC methodology for accurate parameter estimation and prediction, enabling efective measures and infrastructure planning in the region.


Introduction
Extreme rainfall events are no longer a distant possibility but a harsh reality that we must face [1][2][3][4].Tese natural phenomena are becoming more frequent, more intense, and more devastating, afecting millions of people and ecosystems worldwide [3,5].Recent scientifc evidence [1,[3][4][5][6][7][8] claims global warming to be the root cause of these extreme weather events.With rising temperatures and increasing levels of atmospheric carbon dioxide, the hydrological cycle is thrown of balance, leading to more frequent and severe foods and droughts across the globe [1,4,5,9].As per the fndings of the Intergovernmental Panel on Climate Change (IPCC) Report [10], rural, underprivileged communities within developing nations face the highest vulnerability to the impacts of climate change.Te repercussions of climate change involve heightened variability in extreme rainfall patterns, leading to fuctuations in river fows and an increased occurrence of droughts and foods [11][12][13][14].
Like many other countries, Tanzania is vulnerable to these changes [15].For example, according to the Tanzania Climate Status Statement of 2022 [16], on February 26, 2022, extreme rainfall and thunderstorms led to casualties and displacement in Sumbawanga, Nkasi, and Kalambo districts, Rukwa region.Additionally, in Ludewa district (Njombe region), extreme rainfall between February 22 and 23, 2022, caused the destruction of four bridges, disrupting road transport temporarily.Likewise, between April 26 and 28, 2022, extreme rainfall caused severe fooding in the Mbeya and Songwe regions of Southern Highlands Tanzania [16].Te fooding resulted in 5 fatalities, 21 injuries, 5 missing people, and the displacement of 3150 individuals.Numerous buildings, including houses, schools, and religious structures, were brutally damaged.Additionally, agricultural felds spanning 24,700 acres of farm felds, including paddy, maize, and banana farms were destroyed, and a signifcant number of livestock perished [16].Te fooding also led to the destruction of bridges and roads, impacting transportation in various areas [16].Generally, extreme rainfall events have already taken a severe toll on the country's agricultural sector, causing signifcant losses in maize and rice production in the Southern Highlands [15,[17][18][19].Tanzania and other countries that rely heavily on agriculture, the impact of extreme rainfall is already posing great challenges for their economies.Terefore, it is crucial for climatologists, meteorologists, economic planners, and other policymakers to have a better understanding of heavy rainfall patterns and their future behaviors for efcient and efective planning, decision-making, and mitigation purposes [17].
Extreme rainfall is characterized as the highest amount of rainfall experienced within a 24-hour period annually [20].One way to analyze changes in extreme rainfall is to use statistical distribution models.Tese models are useful for studying extreme rainfall and can aid in food mitigation and control [21,22].Numerous probability distributions can be employed to analyze extreme rainfall patterns.However, it is important to note that the data frequently exhibits nonnormal distribution characteristics across multiple regions.While the Anderson-Darling and Kolmogorov-Smirnov tests have often been used for ftting probability distributions, they may not be sensitive enough to identify substantial deviations from an assumed distribution at a particular location [20].To address this issue, researchers suggest utilizing L-moments as a reliable method to estimate higher statistical moments.Te L-moments method allows for reasonable estimates to be obtained even with as small as 20 sample sizes and without making assumptions about the underlying distribution [20,23].Tis method also ofers an excellent ft to the entire cumulative distribution function using only a limited number of parameters.
Several researchers [24][25][26][27][28] have used L-moments to ft probability distributions of extreme rainfall series, providing a powerful alternative for ftting probability distributions of non-normally distributed data, which is particularly important when dealing with extreme rainfall data where signifcant deviations from an assumed distribution can occur at a specifc location.
Tis study aims at fnding the most suitable probability distribution model for the extreme rainfall series utilizing data from Southern Highlands, Tanzania, and to estimate the extreme rainfall amounts that might occur once in 5, 10, 20, 50, and 100 years to come.Te L-moments help guide the selection of candidate distributions and provide insights into the shape of the data.It is used to perform distribution analysis to select the best distributions.Te L-moments method is a popular approach for analyzing the statistical properties of hydrological variables, especially in regions with limited data [29].Tis study applies fve distinct probability distribution models, namely, the GEV, generalized pareto (GPD), generalized logistic (GLO), generalized normal (GNO), and the Pearson type three (PE3) distributions in order to evaluate the best ft model.Tese distributions are chosen due to their simplicity, superior performance, and widespread usage in frequency analysis of extreme events, as highlighted by [20,23,30].It is widely acknowledged that numerous approaches have been extensively employed to estimate the parameters of potential probability distributions, such as the maximum likelihood method, method of moments, L-moments method, and least squares method [22,[31][32][33].Nonetheless, the Bayesian method has not seen widespread utilization.In this research, we utilize the maximum likelihood estimation, L-moments, and Bayesian method to determine the parameters for the chosen probability distribution, and we conduct an extreme value analysis of annual maximum daily rainfall to compare the efcacy of the frequentist approach with the Bayesian approach.Bayesian modeling allows for more supplementary information via prior knowledge, which can improve statistical inference on extremes, considering that extreme data are rare by nature [8].
Overall, it is crucial to model extreme rainfall to understand its potential impacts and inform policy decisions [15,17,34,35].By using appropriate modeling techniques, we can better understand the behavior of extreme rainfall patterns and how often they occur in the Southern Highlands of Tanzania and help mitigate their efects on the economy, agriculture, and society [36].Te rest of this work is structured as follows: Te data, study area, and methodologies for model ftting are described in Section 2 along with the best ftting techniques.Lastly, the statistical analysis results of the most suitable model are presented in Section 3, whereas Discussion and Conclusions are presented in Section 4 and Section 5, respectively.

Study Area.
Figure 1 shows the study area that covers four stations in the Southern Highlands of Tanzania, which is the administrative zone comprised of the four regions: Iringa, Rukwa, Mbeya, and Ruvuma.Te Highlands zone is at latitude 05 °− 11 °S and longitude 30 °− 39 °E.Te rainfall amount in the area varies between 823 mm to 2,850 mm.Te rainy season occurs from November to April.During dry season, higher elevation regions face some light rain and mists from May to August.Most rainfall is attributed to thunderstorms forming over the Lake Nyasa (Lake Malawi), and the areas facing towards the lake tend to receive more rain.Te highlands have lower temperatures ranging from 13 °C to 19 °C than the surrounding lowlands.Te highest elevated regions face nighttime frosts regularly from June to August.

2
Advances in Meteorology As per information provided by the TMA in 2011 and 2019 [15,17], rainfall that exceeds 50 mm within a 24-hour period is classifed as intense (or extreme) and has the potential to result in fooding.Figure 2 presents the original daily rainfall data spanning from 1990-2020, depicting instances of extreme rainfall occurrences exceeding the threshold of 50 mm for each region as indicated by the values above the red line (i.e,rainfall > rbin 50 mm).

Probability Distributions.
Extreme value theory has emerged as a crucial statistical discipline in the feld of meteorology and hydrology, playing a signifcant role over the past six decades.Tis has led to its widespread adoption in numerous research studies.While many probability distributions can be applied to represent empirical rainfall extremes [22], this study focuses on fve commonly used distributions: GPD, GLO, GNO, PE3, and GEV distributions.Te aim is to identify the most suitable distribution that best fts the empirical annual rainfall extremes.Leveraging L-moments, we determine the best ft distribution and utilize it to analyze and predict extreme rainfall patterns in the Southern Highlands Region of Tanzania.By employing maximum likelihood, L-moments, and Bayesian methods to estimate the parameters μ, σ, and ξ, we gained valuable insights into the behavior of extreme rainfall events.Te best ft model facilitates the estimation of return levels for extreme rainfall, which is essential in designing robust infrastructure and efective food protection measures.

Generalized Extreme Value (GEV).
Te GEV distribution arise from a combination of three models, namely, the Gumbel, Frechet, and Weibull whose cumulative distribution function as given by [9,20,22] is as follows: defned on the set of values x: 1 + (ξ(x − μ)/σ) > 0  .Te parameters of this family of models satisfy the conditions − ∞ < μ < ∞, σ > 0 and − ∞ < ξ < ∞.Te distribution in question was independently discovered by both [37,38].Its three parameters, namely, location μ, scale σ, and shape ξ are used to classify it into three types of extreme value distributions, with ξ being the key diferentiating factor.When ξ > 0 and ξ < 0, the distribution falls under Frechet, and Weibull, respectively.When ξ � 0, it corresponds to the limit as ξ ⟶ 0 of equation ( 1), leading to the Gumbel family of distributions.

Generalized Pareto (GPA).
Te cumulative distribution function for GPA is given as [20]: whereby σ, μ, and ξ are scale, location, and shape parameters.Te range of x is given by

Generalized Logistic (GLO).
Te cumulative distribution function for GLO given by [20,39] is as follows: whereby Advances in Meteorology σ, μ, and ξ are scale, location, and shape parameters.Te range of x is given by

Generalized Normal (GNO).
Te probability density function (PDF) and the cumulative probability distribution (CDF) for the GNO distribution are respectively given by where σ, ξ, and μ are the scale, shape, and location parameters, respectively.Te range of x is defned as follows:

Pearson Type Tree (PE3
). PE3 has the following cumulative probability distribution [20]: 2.4.L-Moments Method.L-moments are based on linear combinations of data or statistics that are sorted in ascending order [40,41].Tey are an improvement over normal moments for assessing probability distribution shape and estimating parameters, especially for small sample sizes.According to [42], L-moments can be computed linearly using Shifted Legendre Polynomials and the uniform distribution function as the foundation, providing a more robust estimate for a given amount of data.Te frst four L-moments of a probability distribution can be calculated using a sorted sample vector (X i: n ) of length n, ranked value i ′ th, and expected value function E as follows: Based on the information presented above, it is evident that the frst L-moment represents the mean value of a data set by calculating the expected value of all individual values.Te second L-moment, on the other hand, measures the expected diference between the largest and smallest values in the data set and is equal to half of this expected value.

L-Moments Ratios.
Te second L-moment is always normalized by the mean value, resulting in a L-moment coefcient of variability (L-Cv) calculated as the ratio of the second L-moment to the frst L-moment: Higher L-moment ratios are obtained by scaling the corresponding higher L-moments using the second Lmoment as a reference point [29]: Te L-Cv (τ) is similar to the standard CV and Lskewness (τ 3 ) measures the lack of symmetry in a distribution, while L-kurtosis (τ 4 ) measures the peakedness of the distribution.L-skewness and L-kurtosis are constrained to lie between (− 1,1) and L-kurtosis is bounded by τ 3 to be not less than − 0.25.L-moments can estimate many probability distributions and are computed from probability-weighted moments.More information can be found in [43].

Probability Weighted Moments.
According to [20,39], probability-weighted moments (PWM) are calculated from sorted data values, where x 1: n ≤ . . .≤ x n: n represents the sorted sample of size n.Te PWM of order r is defned by the following equation: where E • { } is the expectation and F(X) is the cumulative distribution function of a random variable X. Te estimators of PWM can be obtained using the equation as given by [20,39]: Te frst four unbiased sample L-moments (estimators) are given by [20] as follows: Te sample L-moment ratios are given by where t and t r are the natural estimators of τ and τ r respectively, and the sample L-Cv is given by

L-Moment Ratio Diagram.
Te L-moment ratio diagram is a helpful tool for analyzing and testing the ft of probability distributions.Tis technique involves plotting L-Skewness (τ 3 ) versus L-Kurtosis (τ 4 ) and is commonly used to choose the most viable probability distribution for a given dataset.In the study by [20], two-parameter distributions were represented as points on the diagram, while threeparameter distributions were depicted as lines.Te best ft distribution was determined by comparing the location of the data point (τ 3 , τ 4 ) to the line representing the candidate distribution and by measuring the distance between (τ 3 , τ 4 ) and (τ 3 , τ Dist 4 ), where (τ 3 ) and (τ 4 ) represent the skewness and kurtosis values of the observed data, and (τ Dist 4 ) is the kurtosis value of the candidate distribution [20].

Block Maxima Technique.
Te extreme values of a set of independent observations, X 1 , X 2 , . . ., X n , were analyzed by splitting the data set into m blocks, each of length n, such that n is a large number.By taking the maximum value Advances in Meteorology within individual blocks, a new set M n,1 , M n,2 , . . ., M n,m is generated, that can be modeled using the appropriate extreme value probability distribution.It is important to choose an appropriate block size, since if n is too small, the limit model may not provide a good approximation, which can lead to biased parameter estimation and unreliable extrapolation.Conversely, when n is extra large, the number of blocked maximum will be insufcient, causing large variances in estimation.Tus, it is essential to fnd a compromise between bias and variance.

L-Moments Estimation Method.
Te L-moments method estimates distribution parameters by comparing the sample L-moments with the corresponding population L-moments [39].L-moments ofer quantifcations of location, dispersion, skewness, kurtosis, and other aspects of probability distributions or datasets.
For more information about L-moments of various probability distribution models, refer to [44].

Maximum Likelihood Estimation Method. MLE aims
to identify model parameters that optimize the likelihood function with respect to those parameters.However, MLE is most efective when sample sizes are large (as n ≥ 30).Despite this, some researchers still prefer to use MLE even with small sample sizes because it can provide accurate estimates for certain values of the shape parameter.[9] commented that MLE is reliable when the shape parameter (ξ) is greater than − 0.5, but its consistency decreases when the shape parameter is less than − 0.5.

Computing Maximum Likelihood Parameter Estimates.
According to [9], if ξ ≠ 0, the log-likelihood function of the GEV with three parameters can be expressed as follows: given that In addition, assuming that X � X 1 , . . ., X n   is a random sample from the GPD with the extreme value X (n) above a threshold u, the log-likelihood function of the GPD can be obtained using the following equation provided by [33,45]: 6 Advances in Meteorology Similary, MLE for other distributions were obtained but not presented in this work.

Bayesian Modeling Method.
In the Bayesian approach, the model parameters may not be random, but they are treated as random by assigning probability distributions to describe the uncertainty of their values [21,46].According to [47], three critical steps are involved in the Bayesian approach to parameter estimation, which include specifying a probability model for the data, constructing the posterior distribution, and checking the model for its outputs before using them for inference.
2.6.1.Fundamental Principles.Let a dataset x � (x 1 , . . ., x n ) representing realizations of a random variable with a density belonging to the parametric family F � f(x; θ): θ ∈ Θ as detailed by [9].Prior beliefs about the parameter θ can be expressed using the probability density function π(θ), independent of the observed data.According to [9], the likelihood for θ is given by: Te combination of prior knowledge and the probability distribution of the data is achieved through Bayes' theorem, leading to the posterior distribution of θ as given by [9] where π(θ) is the prior distribution for θ, L(θ | x) is the likelihood function, and π(θ | x) is the posterior distribution.
In extreme value analysis, the goal is often to predict the probability of extreme events in the future based on observed data.Te predictive distribution, which describes the likelihood of diferent outcomes of a future experiment, can be obtained within the Bayesian framework.If y represents a forthcoming observation with a probability density function denoted as f(y | θ), and π(θ | x) corresponds to the posterior distribution of θ based on the observed data x, the predictive distribution of y given x can be described as follows [9]: Computing this integral may be difcult, but simulationbased techniques like Markov chain Monte Carlo (MCMC) can provide estimates of the posterior distribution through a simulated sample [9].Bayesian procedures can be advantageous when a suitable prior distribution can be specifed [9].Te main objective is to obtain a Markov chain that converges to the posterior distribution in Bayesian statistics.To accomplish this, there are several algorithms available for implementing MCMC simulations, but this study uses the Metropolis-Hastings algorithm.

Return Level Estimates.
After ftting a probability distribution model to the data, the next step is to estimate the return levels for rainfall.Te T-year return level, denoted as x T , is the level that is exceeded on average only once in T years.Te cumulative probability of nonexceedance can be given by the following equation [9,48]: where T is the return period.For instance, in order to obtain x T from the GEV model with return period T, we need to solve for x in equation ( 1) which gives the following [9,48]: In addition, for the GPD model, the return level x m is defned as the extreme level that is exceeded on average once every m observations [9].Tis is given by where ζ u represents the probability that the value X exceeds the threshold u, i.e., ζ u � Pr(X > u).Te quantity x m can also be defned as the m-observations return level [9].

Statistical Analysis.
Te statistical analysis presented in this section reveals some interesting insights about the yearly maximum daily rainfall in four stations.According to the results shown in Table 1, the Ruvuma region receives the highest amount of rainfall, indicating that this region may be more susceptible to fooding and other related hazards.In contrast, Rukwa, the station with the lowest rainfall amount, may experience some challenges related to water scarcity and drought particularly during the dry seasons when water sources dry up.Tese fndings underscore the signifcance of considering local weather patterns and climatic conditions when planning for development and resource management in the region.
Using the block maxima method, Figure 3 displays the annual maxima daily rainfall data collected over a time period of 31 years from four stations in the Highlands Region.Upon closer inspection, the results indicate that there is no clear trend in the data, suggesting that these stations have experienced a relatively stable pattern of extreme rainfall events over the past few decades.

Data Assumptions and Trend Detection Analysis.
Before engaging in extreme value modeling of rainfall data, a thorough examination of the data's underlying assumptions was conducted.Potential trends in the series, capable of infuencing the modeling process, were carefully assessed for normality, homogeneity, and stationarity.Tests such as Shapiro [49,50] and augmented Dickey-Fuller (ADF) [51,52] were employed for this purpose.Given the diverse distribution of stations with varying rainfall patterns and geographical locations, a station-by-station data analysis was conducted, rendering the homogeneity test inapplicable for this study.For more detailed procedures, refer to [11,19,53,54] and [55].Table 2 displays the results, indicating that the data fulfll the mentioned assumptions and can be used for further analysis.Te study also analyzed the Mann-Kendall trend test and Sen's slope estimator for annual maxima rainfall in all stations.Kendall's tau, a measure of the strength and direction of the trend in the data, indicates a weak positive correlation between time and rainfall.All p values are greater than 0.05, suggesting insufcient evidence to reject the null hypothesis of no trend.Sen's slope values are relatively small, indicating a gradual increase in rainfall over time.Te trend analysis results, illustrated in Figure 4, suggest a weak positive trend in rainfall, aligning with the return levels obtained in Table 3. Te fndings are consistent with the study conducted by [56].4, provide valuable insights into the characteristics of extreme rainfall in these regions.Te λ 1 values represent the mean rainfall depth, with Ruvuma having the highest average rainfall depth of 74.3066.Te λ 2 values indicate the variability in rainfall depths, with Iringa having the lowest spread (scale) of 8.34 kurtosis (τ 4 ) and Skewness (τ 3 ) values reveal the heaviness of the tails, and the departure from symmetry, respectively.Iringa shows a slightly right-skewed distribution (0.3177) and a relatively high kurtosis (0.2267).Te λ 5 values refect the departure from symmetry in the upper tail, with Iringa having a thicker upper tail (0.9665).Te τ r values represent the overall shape, and all stations exhibit similar values, indicating moderate peakedness.Tese fndings enable the construction of L-moment ratio diagrams to identify the best-ft distribution for extreme rainfall in each region, facilitating a comprehensive understanding of the rainfall characteristics and aiding in effective hydrological modeling.To ensure the reliability of the results, the L-moments ratio diagram was employed to assess the relative distribution of extreme rainfall.Tis diagram played a crucial role in selecting the most appropriate distribution among GPA, GNO, GEV, GLO, and PE3 across the four stations.Based on the Lmoments ratio diagram results of Figure 5, it was evident that the GEV distribution exhibited the best ft for extreme rainfall in all four stations.Consequently, parameter estimation for the GEV distribution was carried out using MLE, L-moments, and Bayesian MCMC.

GEV Parameter Estimates.
To estimate the GEV distribution, we extracted the blocked maximum of daily annual maximum rainfall data from all four stations, using blocks of n � 365 days to ensure that our samples were sufciently large.We then applied maximum likelihood estimation to ft the model to the N � 31 annual maximum.

Advances in Meteorology
Also, the Bayesian MCMC was used to provide the Bayesian parameters of the annual maxima rainfall using noninformative priors.When selecting noninformative priors for GEV modeling, we minimized prior infuence to allow the data to play a key role in shaping the posterior distribution [57].For the location parameter μ, we deployed a Normal distribution that is broad and symmetric distribution centered at 0. Te scale parameter σ, was modeled with a wide normal distribution to encompass the expected range of extreme values.Similarly, a normal distribution that permits a wide range of shapes, indicating minimal prior information about the tail behavior was used for the shape parameter ξ.Following the methodology outlined by [2,57,58], the joint density of our parameters α, μ, and ξ was represented as Tus for the noninformative priors we employed the following: f α (α) ∼ N(0, 10000), f μ (μ) ∼ N(0, 10000), and f ξ (ξ) ∼ N(0, 10000).Here, N(0, 10000) denotes a normal distribution with a mean of 0 and a variance of 10,000, ensuring a high variance to afrm the absence of external information.
Table 5 shows the outcomes for MLE, L-moments, and Bayesian posterior parameter estimations for the GEV over all four stations.We plotted the diagnostic plots for the GEV model for all three methods as shown in Figures 6-8.Te model diagnostic plots for rainfall at Iringa exhibit a satisfactory ft for the GEV distribution across all methods.Although we conducted similar diagnostic plots for the other stations, we have not included them in this study.Nevertheless, all important diagnostic plots including quantile 10 Advances in Meteorology plot, the probability plot, density plot, and the return level plot present compelling evidence that GEV distribution is an appropriate model that fts the block maxima data in all stations.
From Table 5, the MLE shows that the shape parameter ξ, was estimated to be 0.249 for Iringa, 0.089 for Mbeya, − 0.06 for Rukwa, and 0.0199 for Ruvuma.Te positive values for Iringa and Mbeya indicate that these two regions have an upper bound, while the negative value for Rukwa indicates a lower bound.Te value for Ruvuma is close to zero, indicating that the distribution for this region is approximately unbounded.Notably, the confdence interval for ξ contains 0, indicating that a Gumbel distribution can also provide a better ft for the data.Tis fnding is supported by the profle confdence interval for ξ, which clearly shows that the value ξ � 0 falls well within the interval for all stations.Te values have small diference when L-moments and Bayesian MCMC were used.Te MCMC estimation results revealed the existence of an upper fnite endpoint to the distribution due to positive shape parameter.Tis estimation is in agreement with the MLE approach when used for Iringa station.Te diagnostic plots for GEV estimation indicated in Figures 6-8 shows in all methods the empirical quantiles are linearly related with model quantiles.In both methods simulated data are right-distributed and afected by extreme values.

Trace Plots and Posterior Densities Analysis.
Trace plots and posterior densities obtained from Bayesian MCMC analysis for the GEV distribution provide valuable insights into the parameters of extreme rainfall distribution in the regions.According to Figure 9, the trace plots show the convergence and mixing of the MCMC chains, indicating reliable parameter estimation.Te absence of a trace plot for the shape parameter (ξ), indeed shows that the estimated shape parameter is close to zero.Tis indicates that the data may also be well-ftted by the Gumbel distribution, which is a special case of the GEV distribution with a shape parameter of zero.Te reason behind this behavior is that the shape parameter determines the tail behavior of the distribution.When the estimated shape parameter is close to zero, it suggests that the tails of the data are not signifcantly heavy or light compared to the Gumbel distribution, which has a standard exponential tail.Terefore, the model is efectively converging towards the Gumbel distribution, and the trace plot for the shape parameter may not be displayed.
Te posterior densities of the location parameter (μ), and scale parameter (σ) provide information about the uncertainty associated with these parameters.Based on the trace plots and well-defned posterior densities, we can conclude that the MCMC algorithm has successfully explored the parameter space and converged to stable estimates for μ, σ, and ξ.Tis indicates that the Bayesian MCMC method is efective in estimating the parameters of the GEV distribution for extreme rainfall in Iringa.Similary, fairly good trace plots and posterior densities for the GEV distribution parameters were obtained for other stations but not presented in this work.

Rainfall Return Levels for Diferent Return Periods.
Te return level estimates obtained using diferent methods provide important insights into the potential intensity of extreme rainfall events in Iringa.Te 100-year return level estimated using maximum likelihood estimation (MLE), is projected to be 128.45mm of annual daily maximum rainfall, with a confdence interval ranging from 54.56 mm to 202.344 mm.Tis implies that, without intervention to mitigate climate change, Iringa could experience intense rainfall events of such magnitudes once every hundred years.Te Bayesian MCMC method, which incorporates uncertainty quantifcation, provides a slightly lower estimate for the 100-year return level at 91.83 mm, with a wider confdence interval ranging from 188.40 mm to 584.76 mm.Tis suggests that the Bayesian approach acknowledges the inherent uncertainties in the estimation process and presents a broader range of possible outcomes.Comparatively, the Lmoments method estimates the 100-year return level at 81.20 mm, with a narrower confdence interval ranging from 124.74 mm to 194.38 mm.Tese variations in estimates across methods highlight the importance of considering diferent approaches and their associated uncertainties in assessing extreme rainfall patterns.As the return period increases, both the estimated return levels and the width of the confdence intervals tend to increase, indicating a higher likelihood of more intense rainfall events occurring over longer periods.Table 3 presents the estimations of return levels for diferent return periods in Iringa.Similar return levels for diferent return periods were made for other stations but were not presented in this study.Tese fndings provide valuable information for understanding the potential risks associated with extreme rainfall in the regions and can aid in decision-making for infrastructure planning and climate adaptation strategies.

Performance of MLE, L-Moments, and Bayesian Methods.
Table 6 presents the performance comparison of Bayesian MCMC, MLE, and L-moments, in estimating the parameters of the GEV distribution.Te evaluation is based on two metrics, MAE (mean absolute error) and RMSE (root mean square error).Te results show that the Bayesian MCMC method achieved the lowest MAE (0.312) and RMSE (0.243) values compared to the other two methods.Tis indicates that the Bayesian MCMC method provides estimates that are closer to the true values of the GEV distribution parameters, resulting in smaller prediction errors on average.Te MLE method obtained slightly higher MAE (0.327) and RMSE (0.268) values compared to the Bayesian MCMC method, indicating slightly larger prediction errors.On the other hand, the L-moments method yielded the highest MAE (0.342) but had a lower RMSE (0.253) compared to the MLE  Advances in Meteorology 13 method.Tis suggests that the L-moments method has a higher average prediction error but exhibits less variability in the predictions.Based on these fndings, we can conclude that the Bayesian MCMC method outperforms the MLE and L-moments methods in estimating the parameters of the GEV distribution for the given rainfall data.It provides more accurate and precise parameter estimates, resulting in better predictions of extreme rainfall events.Sample number 6,000 7,000 8,000 9,000 10,000

Discussion and Conclusion
Te statistical analysis of yearly maximum daily rainfall at four stations revealed valuable insights into the rainfall patterns and characteristics in the Highlands Region.Te results demonstrated signifcant diferences in rainfall amounts among the stations, with Ruvuma receiving the highest amount and Rukwa experiencing the lowest.Tis suggests that the Ruvuma region may be more susceptible to fooding and related hazards, while Rukwa may face challenges related to water scarcity and drought during the dry seasons.Te analysis of the annual maxima rainfall data using the Block Maxima method indicated a relatively stable pattern of extreme rainfall events over the past few decades, without a clear trend.Tese fndings imply that the risk of fooding or other related hazards in these regions may remain constant over time.However, it is crucial to continue monitoring weather patterns and extreme events to ensure the implementation of appropriate measures for protecting people and infrastructure from the impacts of climate change.
Te L-moments ratios provided further insights into the characteristics of extreme rainfall.Te values of λ 1 (mean rainfall depth) showed that Ruvuma had the highest average rainfall depth, while Rukwa had the lowest spread (scale) represented by λ 2 .Te kurtosis (τ 4 ) and skewness (τ 3 ) values indicated the heaviness of the tails and the departure from symmetry, respectively.Iringa exhibited a slightly right-skewed distribution and relatively high kurtosis.Te λ 5 values refected the departure from symmetry in the upper tail, with Iringa having a thicker upper tail.Overall, the Lmoments ratios suggested that the extreme rainfall in these regions had moderate peakedness and varying degrees of tail heaviness and asymmetry.
Te L-moments ratio diagrams played a crucial role in identifying the best-ft distribution for annual extreme rainfall in each region.Te analysis showed that the GEV distribution provided the best ft for extreme rainfall in all four stations, as indicated by Figure 5. Te fnding was further supported by the parameter estimation using maximum likelihood estimation (MLE), L-moments, and Bayesian Markov Chain Monte Carlo (MCMC) methods.Tis fnding corroborates previous research by [17,59] highlighting the efectiveness of the GEV distribution in characterizing extreme events.
Te GEV parameter estimates revealed important characteristics of the rainfall distributions in the study area.Te shape parameter (ξ) values indicated whether the distributions had upper or lower bounds or were unbounded.Iringa and Mbeya exhibited positive shape parameters, suggesting an upper bound, while Rukwa had a negative shape parameter, indicating a lower bound.Ruvuma had a shape parameter close to zero, implying an approximately unbounded distribution.Te confdence intervals for ξ contained the value of 0, indicating that a Gumbel distribution could also provide a reasonable ft for the data.Tis fnding is supported by the absence of a trace plot representing the shape parameter (ξ), indicating that the estimated shape parameter is near zero.Tis suggests that the data could be well-suited for ftting the Gumbel distribution, which is a special case of the GEV distribution with a shape parameter of zero.Tis behavior stems from the fact that the shape parameter infuences the tail characteristics of the distribution.When the estimated shape parameter is close to zero, it implies that the data's tails are not signifcantly heavier or lighter compared to the Gumbel distribution, which has a standard exponential tail.
Additionally, the estimation of return levels based on the GEV distribution provided valuable information on the magnitude and frequency of extreme rainfall events in the Highlands Region.Te results of this study provide valuable insights into extreme rainfall patterns and their impacts in the Southern Highlands Region of Tanzania.However, it is crucial to acknowledge that each region may possess unique characteristics that infuence the occurrence and consequences of extreme rainfall.As such, caution should be exercised when generalizing these fndings to other locations.Nonetheless, the methodology and analysis employed in this study can be universally applied as a framework for deducing the best-ft distribution that characterizes extreme rainfall behaviors.By utilizing a rigorous statistical approach and considering relevant climatic factors, this study has contributed to understanding the underlying patterns of extreme rainfall events in Southern Highlands Region of Tanzania.Te fndings can serve as a basis for similar studies in diferent regions, adapting the methodology to the specifc local context.
Furthermore, the identifcation of potential links between extreme rainfall and climate change in Tanzania highlights the need for further research and analysis.To achieve a more comprehensive understanding of the universal aspects of extreme rainfall patterns and their implications for climate change adaptation worldwide, it is essential to undertake broader studies encompassing multiple regions and considering a longer time span.Such studies could explore the infuence of global climate change phenomena on regional extreme rainfall events, facilitating the development of more efective adaptation and mitigation strategies.

Data Availability
Te corresponding author will provide data used in this research upon request.

2. 2 .
Data Description.Rainfall data (in millimeters) was collected from the Tanzania Meteorological Authority (TMA) from four diferent weather stations in the Southern Highlands over a period of 31 years from 1990 to 2020.Each station represents a diferent region in the highlands and has its own unique rainfall patterns and location.Because the data were analyzed station by station, a test for homogeneity was not relevant for this study.

Figure 1 :
Figure 1: Tanzania map showing the Southern Highlands regions.

Figure 6 :Figure 7 :
Figure 6: Diagnostic plots of the MLE method in Iringa station.

Figure
Figure Diagnostic plots of L-moments method in Iringa station.

Figure 9 :
Figure 9: Trace plots and posterior densities of location μ, and scale σ parameters for Iringa station.

Table 1 :
Descriptive statistical analysis of annual maximums.

Table 2 :
Assumptions on extreme rainfall series using diferent tests.

Table 6 :
Performance comparison of MLE, L-moments, and Bayesian MCMC.