Spatial Dependence Analysis of Weekly Moving Cumulative Rainfall for Flood Risk Assessment

Abstract

Climate change has intensified the frequency and severity of extreme weather events, necessitating a nuanced understanding of flood patterns for effective risk management. This study examines flood risk in the Chi watershed, Thailand, using Weekly Moving Cumulative Rainfall (WMCR) data from 1990 to 2021. We employ extreme value copula analysis to assess spatial dependence between meteorological stations in the watershed. Nine bivariate generalized extreme value (BGEV) models were evaluated using the Akaike Information Criterion (AIC) and the Likelihood Ratio test (LRT) to ensure model robustness. The BGEV model revealed higher tail dependence among stations near the bay of the watershed. We also calculated the flood recurrence period to estimate flood events’ frequency and potential severity. Stations ST5 (Khon Kaen), ST6 (Tha Phra Khon Kaen), and ST8 (Maha Sarakham) were identified as potential hotspots, with higher probabilities of experiencing extreme rainfall of approximately 200 (mm.) during the rainy season. These findings provide valuable insights for flood management and mitigation strategies in the Chi watershed and offer a methodological framework adaptable to other regions facing similar challenges.

1. Introduction

The impacts of climate change on river discharges and subsequent flooding have become increasingly significant. Consequently, flood management to minimize the impacts of climate change has become a priority at national and regional levels [1]. Floods are one of the most catastrophic natural events, caused by extreme rainfall that exceeds a particular catchment’s capacity. Extreme weather events are indeed becoming more severe and frequent due to radical climate changes [2,3]. Global climate change is leading to increased frequency and severity of extreme weather events, such as hurricanes, floods, extreme rainfall, and heatwaves [4,5]. These events are causing devastating consequences for human populations and intensifying erosion processes in both mountainous and plain regions. The resulting heavy inundations have significant social and economic impacts on human life and the surrounding environment. In fact, they account for 31% of the world’s economic loss, and flash floods, in particular, are more dangerous as they offer little time for warning the populace and result in more fatalities than typical floods [6,7]. According to the World Atlas of Natural Disaster Risk, Thailand is among the top 10% of countries worldwide at risk of economic loss from flooding. Extreme river flooding is predicted to cause more than 2 million injuries by 2035–2044, and coastal flooding could injure almost 2.4 million persons by 2100 [8,9,10]. Flooding in the Chi watershed, where farming and animal husbandry are critical to the economy, has been a recurrent problem for about two to three times, with significant flooding occurring every two to three years. The most devastating floods occurred in 1978, 1980, 1995, 2000, and 2001, resulting in a profound impact on human life and property. The 2011 flood in Thailand serves as an illustration of how natural disasters disproportionately affect lower-income groups [11,12], and the area was studied 11 floods occurred repeatedly during 2005 and 2021 [13]. Therefore, assessing flood risks in the Chi watershed is critical for flood management and reducing the impact on human life and property.

Despite the inherent danger of living in flood-prone areas, settlements continue to be located along floodplains due to the benefits offered by rivers and their basins drawn to flood-prone regions due to the fertile soil left behind by receding waters, resulting in bountiful crops and economic prosperity [14]. Low-lying flood-prone areas often undergo development for various land uses, driven by their strategic locations and agricultural suitability, despite being less than ideal for development [15]. Numerous researchers have explored different approaches to studying floods. For instance, Dottori [16] produced a new dataset for river flood hazard maps for the European and Mediterranean Basin regions by evaluating present and future river flood risk scenarios to the cost-benefit assessment of different adaptation strategies to reduce flood impacts. Chatzichristaki et al. [17] analyzed the flash flood on Rhodes Island in order to achieve the goals of this research the effective rainfall was estimated by using the Curve Number (CN) method and the flood hydrograph was estimated by using Soil Conservation Service (SCS) synthetic unit hydrograph. On the other hand, fundamental analysis with extreme value theory is another effective method that takes into account the nature of spatial data and can simultaneously describe the extreme value properties [18,19]. However, flood factors are dependent on each other, and spatial dependence of flood events at several locations has been assessed using different areas [20]. The copula function is modeled by the marginal distribution and multivariate dependency, offering flexibility in making the needful adjustment of the marginal and joint probability functions.

Although most researchers primarily focus on analyzing daily rainfall data, the utilization of weekly cumulative rainfall data analysis can offer numerous advantages. For instance, weekly cumulative rainfall data contributes to the assessment of the accuracy and skill of rainfall forecasts, leading to improvements in forecasting methods [21]. As a result, weekly moving cumulative rainfall data emerges as an invaluable resource for applications in water resource management, agriculture, climate change mitigation, and the enhancement of rainfall modeling. Also, the approach facilitates the prediction of weekly cumulative rainfall, enabling informed decision-making in agricultural practices and aiding the selection of crop varieties suited to the prevailing climatic conditions [22]. Therefore, this study aims to assess bivariate flood risk in the Chi watershed by analyzing the maximum weekly moving cumulative rainfall and calculating the joint probabilities of the maximum weekly moving cumulative rainfall using various extreme value copula functions. The performance of these copulas will be evaluated using goodness-of-fit statistical tests, and the best extreme value copula will be selected for further flood risk analysis. The risk indices considered in this study are the return periods for critical thresholds of the maximum weekly moving cumulative rainfall, including primary, joint, and secondary return periods.

The rest of this paper is organized as follows: We start by presenting the central methodology adopted for this study, primarily focusing on the application of extreme value copula in spatial dependence models. Section 2 introduces the geographical region of interest for this research along with the data we have employed. An extensive overview of the materials and techniques utilized for this research is provided in Section 3. In Section 4, we discuss the outcomes of our research, with a particular emphasis on the spatial aspects and the performance of the extreme value copula models. Section 5 is devoted to an in-depth discussion of these findings, while Section 6 offers our conclusions and suggests recommendations based on our research.

2. Study Area and Data

2.1. Study Area

This study focuses on the Chi watershed in Thailand, which is depicted in Figure 1. The watershed is represented by the dark green line and is situated in the northeastern region of Thailand, spanning between the latitudes of ${15}^{\circ }{30}^{\prime }$ N to ${17}^{\circ }{30}^{\prime }$ N and longitudes of ${101}^{\circ }{30}^{\prime }$ E to ${104}^{\circ }{30}^{\prime }$ E (WCPC [23]). The watershed covers a total area of 49,131.92 km${}^2$, with the majority of the area located in 14 provinces.

2.2. Data

This study utilized weekly moving Cumulative Rainfall (WMCR) data obtained from the Thailand Meteorological Department (TMD) using maximum daily rainfall (MDR) from 1990 to 2021 [24]. The maximum WMCR values for all stations ranged from 145.90 to 274.50 mm., while the average WMCR ranged from 18.39 to 27.80 mm. The station ID, name, provincial area, latitude, longitude, and descriptive statistics for each station are provided in Tables S1 and S2 in the Supplementary Material. To provide a clearer picture, the ridgeline plots of WMCR by season in Thailand for five directly affected stations in the Chi Watershed is presented in Figure 2.

3. Materials and Methods

In this section, we describe the methodology used for our analysis, which involved organizing the weekly moving cumulative rainfall (WMCR) data into two separate data frames. We then applied the maximum annual WMCR to either a block maxima model or a generalized extreme value (GEV) distribution. We chose these models because the performance of each bivariate copula extreme analysis varies.

3.1. Marginal Distribution

This section focuses on the marginal distribution used for estimating the maximum rainfall data. While there are several options available, this study primarily focuses on the generalized extreme value (GEV) distribution [25].

The GEV distribution function is as Equation (1):

$${\displaystyle F\left(x\right)=\exp \left(-{\left(1+\xi \left(\frac{x-\mu }{\sigma }\right)\right)}^{-1/\xi }\right),1+\xi \left(\frac{x-\mu }{\sigma }\right)>0,}$$

(1)

where $\mu ,\sigma >0$ and $\xi$ which are location, scale, and shape parameters, respectively.

3.2. Extreme Value Copulas

In this study, we applied copula analysis due to the non-linear correlations present among the weekly maximum cumulative rainfall (WMCR) at different stations. Sklar’s Theorem [26] illustrates how the copula, a function that connects univariate marginal distributions to a multivariate distribution, can be utilized. This analysis hinges on the concept of extreme value copulas, particularly because they yield multivariate extreme value distributions when combined with Unit Frechet margins [20].

We graphically examined the dependence structure and performed spatial correlation analysis using Kendall’s correlation coefficient [27]. From this, we identified the pairs exhibiting the highest spatial correlations for copula function analysis. We then used nine copula distribution functions specified in Table S4 and estimated the parameters using the maximum likelihood method. The model with the lowest LR-test value was selected for representation.

The class of extreme value copulas is particularly important among copula functions. This is because multivariate extreme value distributions can be obtained by using copulas with Unit Frechet margins [20]. In the case of bivariate analysis, the bivariate extreme value distribution denoted by $G_{\ast }(x,y)$ with unit Frechet margins can be expressed as Equation (2):

$$G_{\ast }(x.y)=exp(-{\mu }_{\ast }{[0,(x,y)]}^C),$$

(2)

where ${\mu }_{\ast }{[0,(x,y)]}^C=(\frac{1}{x}+\frac{1}{y})A\left(\frac{x}{x+y}\right)$.

$$A\left(\omega \right)={\int }_0^{1}max(q(1-\omega ),(1-q)\omega )S\left(dp\right),$$

the conditions for measurement of S is that $s{\int }_0^{1}qS\left(dp\right)={\int }_0^1(1-q)S\left(dp\right)=1$.

We also analyze the spatial correlations (spatial dependence) of the data between the ranking data with the Kendall’s correlation coefficient as Equation (3) under the null hypothesis of independence of X and Y [27];

$$\tau =\frac{2}{\widetilde{n(n-1)}}(\sum _{i<j}sgn(x_i-x_j)sgn(y_i-y_j)).$$

(3)

Then, we selected each pair of the highest spatial correlations to analyze the copula function as follows.

$$F(x_1,x_2,\dots ,x_n)=C(F_1\left(x_1\right),F_2\left(x_2\right),\dots ,F_n\left(x_n\right)),$$

when $F_1\left(x_1\right),F_2\left(x_2\right),\dots ,F_n\left(x_n\right)$ are marginal distribution function, then Copula function, C is unique

The specified model as in Table S4 total nine copula distribution functions from bivariate generalized extreme value (BGEV) distribution. Then we estimate the parameters by using the maximum likelihood method (MLE) [28] according to Equation (1), where the parameter estimate is obtained from Equation (4). Finally, the model with the deviance value of the LR-test was presented.

$$\widehat{\theta }=argmaxL\left(\theta \right)=argmax\sum _{i=1}^{n}lnC(u_1,u_2,\dots ,u_n),$$

(4)

when $C(u_1,u_2,\dots ,u_n)=F(F_1^{[-1]}\left(u_1\right),\dots ,F_n^{[-1]}\left(u_n\right))$ is quasi-inverses of $F_1,F_2,\dots ,F_n$ and $u_n\in [0,1]$.

We define $y_k=y_k\left(x_k\right)={\{1+{\xi }_k(x_k-{\mu }_k)/{\sigma }_k\}}^{-1/{\xi }_k}$ for $1+{\xi }_k(x_k-{\mu }_k)/{\sigma }_k>0$ and $k=1,2$, where the marginal parameters are given by $({\mu }_k,{\sigma }_k,{\xi }_k),{\xi }_k>0$. If ${\sigma }_k=0$ then $y_k$ is defined by continuity. In each of the bivariate distributions functions $G(x_1,x_2)$ given below, the univariate margins are generalized extreme value, so that $G\left(x_i\right)=exp(-y_k)$.

Our findings revealed that the extreme value copulas provided a robust framework for understanding the complex interdependence and potential flood risk in the Chi watershed. This work contributes to advancing methods for flood risk analysis, which are critical in the context of climate change and its impacts on extreme weather events.

3.3. Goodness-of-Fit Statistical Tests

In order to assess the performance of the marginal and joint probabilities, hypothesis tests are used to determine the goodness-of-fit of the collected data to a particular distribution. Specifically, the Kolmogorov-Smirnov are commonly used for univariate analysis and Likelihood-ratio tests or deviance statistic (D) while the bivariate case in this study employs these tests to evaluate the performance of the joint probabilities [25].

3.4. Selection of Extreme Value Copula Models

This study selects the best EVC models by Akaike Information Criterion (AIC) is to select the model that minimises the negative likelihood penalised by the number of parameters as specified in the equation $AIC=-2logL\left(\widehat{\theta }\right)+2k$ and Bayesian information criterion (BIC) for candidate model is defined as; $BIC=-2lnL({\widehat{\theta }}_k|y)+kln\left(n\right)$ where $\theta =$ the set(vector) of model parameter; $L\left(\widehat{\theta }\right)=$ the likelihood of the candidate model given the data when evaluated at the maximum likelihood estimate of $\theta$; k is number of estimated parameters in the candidate model.

3.5. Risk Analysis

Flood risk is the potential harm arising from the interaction between flood-inducing factors and the environment. It’s often represented by the flood return period, which indicates the frequency of a flood of a certain magnitude. The longer the return period, the more severe the potential flood damage. The flood recurrence period is calculated as $T=\frac{1}{1-P}$, where P is the probability that an event won’t exceed a set threshold. The trend of return period T aligns with the probability P. The higher the value of P, the longer the return period, implying a higher potential loss and thus greater flood risk.

The paper suggests that the joint (primary) return periods can be further characterized by “AND” cases. This implies that the joint return period is calculated by considering the probability of both events occurring simultaneously [29,30]. For instance, the joint return period for a flood event could be calculated by considering the probability of both the rainfall event and the river discharge event exceeding certain thresholds within a given time frame as Equation (5):

$$T_{u_1,u_2}^{AND}={\displaystyle \frac{1}{1-u_1-u_2+C_{U_1{U}_2}(u_1,u_2)}},$$

(5)

where $\mu$ indicates the mean inter-arrival time of the two consecutive extreme events. $u_1$ and $u_2$ are specific values of $U_1$ and $U_2$, respectively;$U_1=F_1\left(x_1\right)$ and $U_2=F_2\left(x_2\right)$; $F_1$ and $F_2$ are the cumulative distribution function (CDF) which the Figure 3 shows clearly about the risk area concept. In addition, small circles are scatter points of rainfall(mm.) between ST5 and ST6, dash lines are threshold for 95% quantile rainfall(mm.) for both stations are shown in Figure 3.

The intensity of flood risk is usually determined by the period of future re-occurrence of heavy inundation, commonly known as the return period. Under consideration of the service time (i.e., n) of a hydrological infrastructure, the risk of failure associated with the return period of a flooding event can be expressed as Equation (6) [30]:

$$R=1-{(1-p)}^n=1-q^n=1-{\left(1-{\displaystyle \frac{1}{T}}\right)}^n,$$

(6)

where R is the risk of failure, p and q are the probability of an event exceeding and not exceeding the set threshold, and T is the return period. The failure probability (R) is the measure of flood risk ranging between 0 and 1, with a greater value indicating a greater risk magnitude. A risk framework that considers more than one variable may provide more support for actual flood control than the conventional analysis [31]. It is essential to characterize floods through multiple variable aspects by a process called bivariate risk analysis. Bivariate risk analysis plays a significant role in taking non-structural safety measures and developing flood mitigation strategies. Bivariate risk analysis in relation to joint return period in AND is defined as Equation (7):

$$R_{u_1,u_2}=1-{\left(1-{\displaystyle \frac{1}{T_{u_1,u_2}^{AND}}}\right)}^n.$$

(7)

In addition, the flow chart of this study are shown in Figure 4 as below:

4. Results

4.1. Marginal Probability Distribution

The parameter estimation using the maximum likelihood estimator (MLE) is used. It was found that the data were suitable for the GEV distribution, as observed from the p-value of KS test and AD test greater than 0.05, which from Table S3, the p-value were between 0.7–0.9. The range of location parameter is (60.21, 99.99), scale parameter is (13.61, 26.55) and shape parameter is (−0.43, 0.33).

4.2. Extreme Copula Value Analysis

In this section Bivariate Generalized Extreme Value (BGEV) distribution is used to model the tail dependence between stations which the component wise block maxima approach is considered for BGEV distribution.

4.2.1. Correlation between Stations

The Annual Maximum Weekly Moving Cumulative Rainfall (AMWMCR) data are presented in the following Figure 5. The Kendall tau correlation is also estimated and compared with estimates of extreme dependence. For convenience, the estimated Kendall tau is presented in the corresponding sub-sections and shows 5 stations that are directly affected in the Chi Watershed, Thailand. Approaches suggested to interpreting the correlation coefficient range is between −1 to 1. We provide three categories; (1) 0.01 to 0.39 is a weak positive correlation,( 2) 0.40 to 0.69 is a medium positive correlation and (3) 0.7 to 1.00 is a strong positive correlation [32]. On the contrary, (1) −0.01 to −0.39 is a weak negative correlation, (2) −0.40 to −0.69 is a medium negative correlation and (3) −0.7 to −1.00 is a strong negative correlation.

The correlation between AMWMCR of two margins (stations) estimated by Kendall tau. Figure 5 right hand side showed the value of correlation among which were −0.2 to −0.4 (low to mild correlation) and the distribution of AMWMCR of all station is right heavy-tailed distribution.

4.2.2. Copula Function Fitting

We estimated some upper quantile probability from fitted BGEV models. The threshold are fixed at $95\%$ quantile of WMCR data, because the value of $90\%,95\%$ and $99\%$ quantile are high enough by the mean residual life plot. In particularly, these amounts are similar to the heavy rainfall criteria from the TMD. All models were considered to analyse the dependence and strength estimated between stations. The summary of dependence estimates were presented as

Table 1. Extreme value copula model for the best BGEV model.

Station ID		$\mathit{\tau }$	Copula Function	D of LR-Test	AIC
ST5	ST6	0.497	Negbilog	0.090	603.695
	ST8	−0.002	Log	0.060	609.535
	ST10	0.062	HR	0.044	599.181
	ST11	−0.058	Alog	0.086	609.767
ST6	ST8	0.060	Log	0.005	630.730
	ST10	0.132	Bilog	0.974	618.780
	ST11	−0.046	Neglog	0.654	634.101
ST8	ST10	−0.221	Neglog	0.001	615.923
	ST11	0.262	Alog	1.177	623.090
ST10	ST11	−0.277	Neglog	0.029	619.134

.

The largest correlation estimated by Kendall tau, between stations ID ST5 and ST6, is 0.497. The strength of dependence is estimated by negbilog copula function. For other stations mostly the strength of dependence is estimated by alog, amix and negbilog, respectively as in Table 1 and Figure 6.

4.3. Risk Analysis

The utilization of multiple flood variables in flood risk analysis provides valuable information for a deeper comprehension of flood characteristics. To reflect historical flooding, return periods were calculated using the BGEV model (Table 1). The risk value of rainfall and joint return periods were presented in

Table 2. Joint return periods from the joint distribution and flood risk for univariate and bivariate distributions.

Region	RP (Years)	Rainfall by Region (mm.)		${\mathit{T}}_{\mathit{and}}$ (Years)	Risk Value
		Region 1	Region 2
ST5–ST6	2	83.80	90.90	3.52	1.00
	5	110.24	117.42	23.57	0.75
	10	122.16	120.93	56.08	0.44
	25	148.67	141.10	211.81	0.14
	50	156.60	147.90	431.69	0.07
	100	160.35	153.45	712.76	0.04
ST5–ST8	2	84.80	74.70	3.50	1.00
	5	109.68	94.38	11.56	0.94
	10	134.73	118.63	44.65	0.52
	25	152.38	139.89	118.12	0.24
	50	172.19	144.39	205.60	0.14
	100	186.30	146.35	294.60	0.10
ST5–ST10	2	84.80	97.30	4.26	1.00
	5	109.68	114.32	12.18	0.94
	10	134.73	127.02	30.54	0.66
	25	152.38	193.00	164.45	0.18
	50	172.19	196.15	240.45	0.12
	100	186.30	198.73	314.59	0.10
ST5–ST11	2	84.80	88.35	3.48	1.00
	5	109.68	119.68	13.84	0.91
	10	134.73	148.85	41.32	0.54
	25	152.38	168.50	77.52	0.34
	50	172.19	172.20	115.09	0.24
	100	186.29	172.20	147.98	0.20

and Figure 7 (see Tables S4 and S5 in the supplementary material for further details). The risk of rainfall was computed using Equation (7) and demonstrated a strong correlation between station and rainfall risk value, as evident in Tables S6 and S7 in the supplementary material.

We aim to analyze flood risk based on extreme rainfall events between two meteorological stations. We focus on extreme rainfall scenarios, the univariate return period (RP) levels are taken to be 2, 5, 10, 25, 50, and 100 Years, the corresponding values of Rainfall by region are calculated, respectively, and the corresponding two-dimensional Copula function values are calculated by the optimal Copula functions of different characteristic variables. According to Equation (5), the corresponding bivariate joint return periods at a given univariate return period level are calculated as flood risk for regions. The computed values are shown in Table 2. The probability distributions of rainfall by region are plotted and corresponding recurring contours are added to the measured data for comparison, as shown in Figure 8 and Figure 9.

Table 2, shows the bivariate joint return period ($T_{and}$ (Years)), and the risk value which is the probability of extreme rainfall in each pair of meteorological stations. The results indicate extreme rainfall by region analyzed the risk of flooding over a range of return periods, including 2, 5, 10, 25, 50, and 100 years, and plotted the corresponding probabilities of exceedance, as shown in Figure 7. Taking the region ST5–ST11 as an example, in the 100-year return periods the estimate of rainfall in ST5 and ST11 are 186.29 mm. and 172.20 mm., respectively on the joint return period 147.98 years making the flood risk of ST5–ST11 by the probability 0.20.

5. Discussion

In this investigation, we implemented copula functions to encapsulate the interdependence structure between the weekly moving cumulative rainfall (WMCR) from two separate stations within the Chi watershed. Our focus was primarily on extreme value copulas, as these have a notable significance in modeling multivariate extreme value distributions. To gauge the performance of marginal and joint probabilities, we executed goodness-of-fit tests. The most suitable copula function was then determined according to the results of these tests. The outcomes demonstrated that the stations directly influenced within the Chi watershed exhibit high-risk values and correlation. In addition, it was observed that there is an escalated risk of extreme rainfall during the rainy season, notably at stations ST5 (Khon Kaen), ST6 (Tha Phra Khon Kaen), and ST8 (Maha Sarakham), which have a greater likelihood of experiencing extreme rainfall in the vicinity of 200 mm. In summary, our findings offer vital insights into flood risk within the Chi watershed. This knowledge is instrumental in informing and enhancing flood management practices and disaster response strategies.

6. Conclusions

This study’s findings can be distilled into three main points. First, the BGEV model demonstrated satisfactory performance in estimating the tail dependence between meteorological stations. Stations located near the watershed’s bay exhibited strong tail dependence, while those geographically distant showed nearly no dependence. Moreover, the BGEVD model’s estimated upper quantile probability exceeded the empirical probability, suggesting its superior capacity in predicting heavy rainfall risk. Second, the research revealed that flood risk analysis can be characterized by extreme rainfall across two meteorological stations. We considered four extreme rainfall patterns 120, 150, 180, and 200 mm. for 2, 10, and 25 weekly events. During the rainy season, stations ST5 (Khon Kaen), ST6 (Tha Phra Khon Kaen), and ST8 (Maha Sarakham) warrant close monitoring due to their heightened probability of experiencing extreme rainfall of around 200 mm. Lastly, the methodology employed in this research can be adapted for analyzing numerical model outputs of climate systems, enabling a comparison with a model fitted to climate observations. Future research should consider significant covariates that impact temperature, such as topographical aspects and proximity to the coast. Modeling a climate variable’s complete distribution for a spatial field may prove more accurate and valuable than focusing solely on extremes. A time-dependent MSP model incorporating climatic covariates could be beneficial for this purpose. The challenges and potentials inherent in this area necessitate further collaborative exploration between climate scientists and statisticians.