Development of extreme value models for precipitation forecasting using a nonstationary process for critical environment sustainable management

1. Introduction

Greenhouse gases in the Earth"s atmosphere have increased rapidly, especially since the Industrial Revolution. Over the past five decades, human activities have caused the concentration of carbon dioxide to surge from 280 ppm, a level that had remained stable for millions of years. In 2022, the concentration reached 421 ppm, which is 50% higher than pre-industrial levels. This has intensified global warming, leading to changes in the climate from its original state and triggering cascading effects on ecosystems and livelihoods, potentially causing a variety of impacts. Thailand is inevitably affected by global warming, and the impacts of "climate change" are becoming increasingly evident. The recurring severe flooding appears to be worsening over time. Statistics from the Meteorological Department of Thailand over the past 10 years indicate that more severe thunderstorms have affected many parts of the country.

Over the past 30 years, Thai communities have frequently experienced natural disasters. The severity of these natural disasters has destroyed houses and damaged farms. Heavy precipitation can cause floods and landslides, particularly in flood-prone areas. Precipitation forecasting is crucial to prepare for, prevent, and address the potential damage from such events. The forecast results provide information about future precipitation, which is beneficial for resource planning to enable effective preventive measures. Many researchers have developed models to forecast precipitation. For instance, Barman et al. (2021) applied Seasonal Autoregressive Integrated Moving Average (SARIMA) and Autoregressive Integrated Moving Average (ARIMA) models to analyze the monthly time series of precipitation in Assam from 1901 to 2017. In 2022, Zhao et al. studied an ARIMA model combined with a Radial Basis Function (RBF) network model to forecast monthly precipitation in Nanchang, Jiangxi Province. Zhao et al. (2023) proposed models combining Ensemble Empirical Mode Decomposition (EEMD), a Long Short-Term Memory (LSTM) neural network, and ARIMA for monthly precipitation forecasting using data from Luoyang spanning 1973 to 2021. Bora and Hazarika (2023) studied the forecasting patterns of precipitation distribution and proposed an ARIMA model to forecast precipitation in Assam and Meghalaya, India. Khan et al. (2023) used an ARIMA model to forecast precipitation in the Klang River Basin, Selangor. Chandran et al. (2023) studied 34 years of precipitation data (from 1976 to 2009) from 10 sub-basins of the Vaigai River in Tamil Nadu using an ARIMA model to identify the appropriate model for precipitation forecasting. Abd-Elhamid et al. (2024) used precipitation data from satellite and ground-based photographs and proposed an ARIMA model for forecasting precipitation in Syria. Sharma and Singh (2024) built an ARIMA model to forecast precipitation over Mandi district, Himachal Pradesh, India. Hendri et al. (2024) studied precipitation data exhibiting seasonal variations using two Holt-Winters models: the additive model and the multiplicative model.

The aforementioned studies focused on developing a model to forecast short-term precipitation. However, spatially accurate climate forecasting, particularly for the long term, remains both crucial and challenging, since it enables relevant agencies to plan effective long-term responses. Many researchers have applied extreme value theory (EVT) to develop models for forecasting long-term extreme values. Analyzing this data with extreme values is crucial for preparing for large-scale catastrophic events, helping to reduce potential losses if such events occur in the future. Therefore, accurate analysis of extreme value data is essential for obtaining precise results and for providing preliminary information for further management. Many researchers have used EVT to develop models for forecasting long-term extreme precipitation, including Chutiman et al. (2019), Busababodhin et al. (2021), Guayjarernpanishk et al. (2023), Sopipan (2023), Chaubey et al. (2023), Lymperi and Varouchakis (2024), and Acero et al. (2024). These studies used data collected from weather stations for analysis. Limitations in data acquisition from these stations may arise due to factors such as the distance between stations and data loss during collection. The data analysis methods in the aforementioned studies apply EVT under a stationary process.

A review of existing literature highlights critical gaps in precipitation forecasting. While many studies focus on short-term predictions or apply EVT for long-term extremes, they often assume a stationary process, which is inconsistent with inherently non-stationary climate data, and primarily rely on weather station data, which suffers from spatial sparsity and potential data loss. This research aims to bridge these gaps by developing a long-term extreme precipitation forecasting model that explicitly incorporates non-stationary processes within an EVT framework and leverages satellite data to overcome ground-based measurement limitations, thereby providing more accurate, spatially refined, and reliable information for climate change adaptation and disaster risk management in Thailand.

2. Materials and Methods

In this research, parameter estimation and long-term precipitation forecasting by employing nonstationary process modeling combined with EVT were investigated. This section presents details of the study area, data preprocessing steps, and the methodology, such as non-stationary process and return level.

2.1 Study area

Thailand is located between latitudes 5°37′N and 20°27′N, and longitudes 97°22′E and 105°37′E. The country is divided into six regions: Northern region (9 provinces, marked in gray), Central region (22 provinces, light brown), Northeastern region (20 provinces, red), Eastern region (7 provinces, yellow), Western region (5 provinces, brown), and Southern region (14 provinces, pink).

Northernmost point: Latitude 20°27"30”N, Mae Sai District, Chiang Rai Province.

Southernmost point: Latitude 5°37′N, Betong District, Yala Province.

Easternmost point: Longitude 105°37"30”E, Phibun Mangsahan District, Ubon Ratchathani Province.

Westernmost point: Longitude 97°22"E, Mae Lan Noi District, Mae Hong Son Province (Fig. 1).

Close

Fig. 1.

Study area presenting the six regions of Thailand.

Export to PPT

2.3 Methods

2.3.1 Generalized extreme value distribution

The Generalized Extreme Value (GEV) distribution, first introduced by A.F. Jenkinson in 1955, was used for analyzing extreme values. This method is appropriate for data collected over regular time intervals, such as weekly, monthly, quarterly, or yearly periods. The data selection process involves selecting the highest value within each time period, based on the following concept:

Let X_i (i = 1, 2, …, n) represent independent variables with the same distribution function, $F(x)$ and define X_(n) = max (X₁, X₂, …, X_n). In the GEV distribution, three parameters are considered: $\mu$ (Location parameter), $\sigma$ (Scale parameter), and $\xi$ (Shape parameter).

In 1978, Galambos developed the Cumulative Distribution Function (CDF) of the GEV distribution for $-\infty <x<\infty$ as follows:

(1)

$G(x;\mu ,\sigma ,\xi )=\exp \left\{-{\left[1+\xi \left(\frac{x-\mu }{\sigma }\right)\right]}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\xi $}\right.}\right\}$

As defined on $\left\{1+\xi \left(\frac{x-\mu }{\sigma }\right)>0\right\}$ when $-\infty <\mu ,\xi <\infty$ and $\sigma >0$

If $\xi =0$ or $\xi \to 0$, we obtain

(2)

$G(x;\mu ,\sigma ,\xi )=\exp \left\{-\exp \left[-\left(\frac{x-\mu }{\sigma }\right)\right]\right\}$

From eq. (1), if $\xi <0$ , the GEV distribution is called the Weibull Distribution, and if $\xi >0$, it is called the Fréchet Distribution. From eq. (2), when $\xi \to 0$, it is called a Gumbel Distribution. The probability distribution function (PDF) of the GEV distribution is as follows:

(3)

$\begin{array}{l}g(x;\mu ,\sigma ,\xi )\\ =\left\{\begin{array}{c}\frac{1}{\sigma }{\left[1+\xi \left(\frac{x-\mu }{\sigma }\right)\right]}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\xi $}\right.-1}\exp \left\{-{\left[1+\xi \left(\frac{x-\mu }{\sigma }\right)\right]}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\xi $}\right.}\right\};\xi \ne 0\\ \frac{1}{\sigma }\exp {\left(\frac{x-\mu }{\sigma }\right)}^{-1}\exp \left\{-\exp \left[-\left(\frac{x-\mu }{\sigma }\right)\right]\right\};\xi =0\end{array}\right.\end{array}$

2.3.2 The GEV distribution under a non-stationary process

Most meteorological, insurance, or economic data analyzed for extreme values are influenced by other variables, such as changes over time or underlying trends, making the data use a non-stationary process. This type of analysis differs from the analysis of data under a stationary process. The characteristics of data under a non-stationary process have been shown in Fig. 2.

Close

Fig. 2.

Maximum daily exchange rate of Euro to Thai Baht from 2009 to 2023.

Export to PPT

The GEV distribution model under a non-stationary process, used to describe the distribution of X_t at time t (where t = 1, 2, .., m), is as follows (Guayjarernpanishk et al., 2021): $X_t\sim \text{GEV}\left(\mu (t),\sigma (t),\xi (t)\right)$

The possible models of non-stationary process are as follows:

Model 1 : $\mu ,\sigma$ and $\xi$ are constants.

Model 2 : $\mu (t)={\beta }_0+{\beta }_{1}t$, where $\sigma$ and $\xi$ are constants.

Model 3 : $\mu (t)={\beta }_0+{\beta }_{1}t+{\beta }_{2}t{}^2$, where $\sigma$ and $\xi$ are constants.

Model 4 : $\mu (t)={\beta }_0+{\beta }_1\exp (-{\beta }_{2}t)$, where $\sigma$ and $\xi$ are constants.

Model 5 : $\sigma (t)=\exp ({\beta }_0+{\beta }_{1}t)$, where $\mu$ and $\xi$ are constants.

Model 6 : $\mu (t)={\beta }_0+{\beta }_{1}t$ and $\sigma (t)=\exp ({\beta }_0+{\beta }_{1}t)$, where $\xi$ is a constant.

Model 7 : $\mu (t)={\beta }_0+{\beta }_{1}t+{\beta }_{2}t{}^2$ and $\sigma (t)=\exp ({\beta }_0+{\beta }_{1}t)$, where $\xi$ is a constant.

Model 8 : $\mu (t)={\beta }_0+{\beta }_1\exp (-{\beta }_{2}t)$ and $\sigma (t)=\exp ({\beta }_0+{\beta }_{1}t)$, where $\xi$ is a constant.

If the result of the parameter analysis follows Model 1, the data used for analysis can be considered stationary. However, if the result follows other models, the data can be classified as non-stationary.

2.3.3 Parameter estimation of the GEV distribution using a non-stationary process

The steps for estimating the GEV distribution parameters using the maximum likelihood estimation (MLE) method are as follows (Guayjarernpanishk et al., 2021):

1. Construct the likelihood function of the probability distribution function of the GEV distribution, which results in $L(\text{β})={\displaystyle \prod _{t=1}^{m}g(x_t;\mu (t),\sigma (t),\xi (t))}$,

where $\text{β}$ is a vector of parameters ${\text{β}}_i$ and $g(x_t;\mu (t),\sigma (t),\xi (t))$ is the probability distribution function of the GEV distribution, with $\mu (t),$$\sigma (t)$ and $\xi (t)$ representing the parameters at $x_t$.

• 2)

  Construct the log-likelihood function of the probability distribution function of the GEV distribution obtained in step 1 and equation (3) for t = 1, 2, ..., m as follows:

• 2.1)

  In the case when $\xi \ne 0$

  (4)

  $l(\text{β})=-{\displaystyle \sum _{t=1}^m\left\{\begin{array}{l}\log \sigma (t)+\left(1+\frac{1}{\xi (t)}\right)\log \left[1+\xi (t)\left(\frac{x_t-\mu (t)}{\sigma (t)}\right)\right]\\ +{\left[1+\xi (t)\left(\frac{x_t-\mu (t)}{\sigma (t)}\right)\right]}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\xi (t)$}\right.}\end{array}\right\}}$

Defined on $1+\xi (t)\left(\frac{x_t-\mu (t)}{\sigma (t)}\right)>0$

• 2.2)

  In the case when $\xi =0$

  (5)

  $l(\text{β})=-{\displaystyle \sum _{t=1}^m\left\{\log \sigma (t)+\left(\frac{x_t-\mu (t)}{\sigma (t)}\right)+\exp \left[-\left(\frac{x-\mu }{\sigma }\right)\right]\right\}}$

Examples of the GEV distribution under a non-stationary process:

When $\mu (t)={\beta }_0+{\beta }_1\exp (-{\beta }_{2}t)$,$\sigma (t)$ and $\xi (t)$ are constants, we obtain: $\begin{array}{l}l\left(\mu (t),\sigma ,\xi \right)\\ =-{\displaystyle \sum _{t=1}^m\left\{\begin{array}{l}\log \sigma +\left(1+\frac{1}{\xi }\right)\log \left[1+\xi \left(\frac{x_t-({\beta }_0+{\beta }_1{\mathrm{e}}^{-{\beta }_{2}t})}{\sigma }\right)\right]\\ +{\left[1+\xi \left(\frac{x_t-({\beta }_0+{\beta }_1\exp (-{\beta }_{2}t))}{\sigma }\right)\right]}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\xi $}\right.}\end{array}\right\}}\end{array}$

• 3)

  Estimate the parameters $\mu (t)$, $\sigma (t)$, and $\xi (t)$, $(\widehat{\mu }(t),\widehat{\sigma }(t),\widehat{\xi }(t))$), by calculating a partial derivative of the functions obtained in step 2, as follows: $\begin{array}{c}\frac{\partial l}{\partial \mu (t)}\left(\mu (t),\sigma (t),\xi (t)\right)=0,\\ \frac{\partial l}{\partial \sigma (t)}\left(\mu (t),\sigma (t),\xi (t)\right)=0\\ \text{and}\frac{\partial l}{\partial \xi (t)}\left(\mu (t),\sigma (t),\xi (t)\right)=0\end{array}$

2.3.4 Model selection

The purpose of model selection is to identify the most appropriate model for the data being analyzed. Assessing model fit is crucial when multiple models result from the analysis. When the data follows a single stationary process, comparing different models to identify the most appropriate one is commonly done using Diagnostic Plots, including the Quantile Plot, Probability Plot, Density Plot, and Return Level Plot. Alternatively, the Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC) can be used.

When selecting an appropriate model for data under a stationary or non-stationary process, using the MLE method generally results in multiple related models being obtained. A model with this characteristic is referred to as a Nested Model. The statistic used for testing in this context is the Deviance Statistic, denoted as D (Guayjarernpanishk et al., 2021).

Define $M_0$ and $M_1$ as the initial model and the model to be compared, respectively, with the condition that $M_0\subset M_1$. The following Hypothesis testing:

H₀: The initial model is appropriate.

H₁: The model to be compared is appropriate.

Using the statistic D, D is defined as follows:

(6)

$D(0,1)=2\{l_1(M_1)-l_0(M_0)\}$

where $l_0(M_0)$ and $l_1(M_1)$ are the maximum log-likelihood values ​​of the models $M_0$ and $M_1$, respectively.

From equation (6), we can see that the statistic D converges in distribution to a chi-square distribution with degrees of freedom equal to k $({\chi }_k^2)$, where k is the difference in the number of parameters between models $M_1$ and $M_0$. Model $M_0$ will be rejected at the level of significance ($\alpha$) if D > $c_{\alpha }$, where $c_{\alpha }$ is the quantile at $(1-\alpha )$ of ${\chi }_k^2$. Therefore, if the model $M_0$ is rejected, it implies that the model $M_1$ better explains the variability in the data than the model $M_0$. The test using statistic D is performed one pair at a time. Once an appropriate model is identified, it is tested against other models until the most appropriate one is found.

Therefore, based on the possible models of non-stationary processes above, when $M_i$ represents the i-th model, we obtain $M_1\subset M_3\subset M_7$, $M_1\subset M_2\subset M_6$,$M_1\subset M_4\subset M_8$ and $M_1\subset M_5$. The sequential pairwise comparison steps are shown in Fig. 3. Additionally, the flowchart outlining the steps for pairwise comparison in the case $M_1\subset M_2\subset M_6$ has been presented in Fig. 4.

Close

Fig. 3.

The sequential pairwise comparison steps.

Export to PPT

Close

Fig. 4.

The steps for pairwise comparison in the case $M_1\subset M_2\subset M_6$.

Export to PPT

Once the most appropriate model is obtained, it will be used for return level estimation at time T years as follows (Guayjarernpanishk et al., 2021):

(7)

${\widehat{R}}_T=\left\{\begin{array}{c}\widehat{\mu }(t)-\frac{\widehat{\sigma }(t)}{\widehat{\xi }(t)}\left\{1-{\left[-\log \left(1-\frac{1}{T}\right)\right]}^{-\widehat{\xi }(t)}\right\};\xi \ne 0\\ \widehat{\mu }(t)-\widehat{\sigma }(t)\log \left\{-\log \left(1-\frac{1}{T}\right)\right\};\xi =0\end{array}\right.$

3. Results

The data used in this research consisted of the daily maximum accumulated precipitation from 208 points, between January 2009 and March 2024. Descriptive statistics of the top two daily maximum accumulated precipitation in each region of Thailand have been shown in Table 1.

As shown in Table 1, the points with the highest accumulated precipitation over the past 15 years, exceeding 100 mm, were located in the Central, Western, North-eastern, and Southern regions. The highest accumulated precipitation was recorded in Nakhon Si Thammarat Province, Southern region, in January 2012, while the second highest was in Tak Province, Western region, in August 2019.

In statistical modeling using EVT, the GEV distribution was applied by setting a monthly block time and analyzing the daily maximum accumulated precipitation data for each month. This analysis utilized the method for calculating extreme values under a non-stationary process across eight models, with data from 208 grid points. The results showed that Model 2 was the appropriate model for three grid points, while Model 1 was the appropriate model for the remaining 205 grid points. The estimated parameters for the appropriate model at the grid points with the top two daily maximum accumulated precipitation values in each region have been presented in Table 2, and the estimated parameter values for the appropriate model for Model 2 at all 3 grid points have been shown in Table 2.

Once the appropriate model for each grid point was obtained, it was used to estimate return levels for time periods of 2, 5, 10, 25, and 50 years. The return levels of daily maximum precipitation in areas that previously experienced the highest precipitation in each region, according to Model 1 in Table 2, and those for Model 2 in Table 3.

Tables 4 and 5 present the return levels of daily maximum precipitation over a period of T years with a probability of 1/T. For example, a return level of 2 years indicates that within a 2-year period, daily maximum precipitation at the forecasted level is expected to occur at least once, with a probability of 1/2 or 0.5. The forecasted values from all 208 grid points for return levels of 2, 5, 10, 25, and 50 years were used to create contour graphs through GIS Kriging interpolation, with the results presented in Fig. 5.

Close

Fig. 5.

Estimated values for the return levels of daily maximum precipitation in the following years: (a) 2 years; (b) 5 years; (c) 10 years; (d) 25 years, and (e) 50 years.

Export to PPT

Fig. 5 indicates that the southern region of Thailand experiences higher daily maximum precipitation compared to other regions. Provinces along the eastern coast, adjacent to the Gulf of Thailand, receive higher daily maximum precipitation than those along the western coast, adjacent to the Andaman Sea. The northeastern region bordering Laos records higher daily maximum precipitation than other parts of the country.

4. Discussion

After statistical modeling with EVT using the GEV distribution under a non-stationary process across eight models, Model 2 was found to be the appropriate model for 3 grid points, while Model 1 was the appropriate model for the remaining 205 grid points. The forecast of daily maximum precipitation recurrence also indicates that the southern region of Thailand is expected to experience higher daily maximum precipitation compared to other regions. This can be attributed to the influence of monsoons from three different directions. From October to January, the southern region is influenced by the Northeast Monsoon, which carries water vapor and moisture from the South China Sea and the Gulf of Thailand to the mainland. This causes heavy precipitation along the eastern coast before it reaches the western coast. In February, the Southeast Monsoon begins to carry water vapor, humidity, and heat from the equator through the Gulf of Thailand to the mainland. This monsoon affects the southern region from February to April, coinciding with the summer season, resulting in high temperatures and precipitation. When the rainy season arrives, the Southeast Monsoon shifts direction to become the Southwest Monsoon, which brings water vapor and moisture from the Indian Ocean and the Andaman Sea to the mainland. This monsoon affects the southern region from May to September, causing heavy precipitation to begin on the western coast before spreading to the eastern coast. Due to the influence of monsoons from all three directions, the southern region experiences heavy precipitation nearly year-round.

In the northeastern region, areas bordering Laos receive more precipitation than other parts of the region. This is due to the rainy season, which lasts from May to October each year and is influenced by the Southwest Monsoon. However, the Dong Phayayen and San Kamphaeng mountain ranges block the winds and precipitation from the southwest monsoon, resulting in lower precipitation in areas located behind the mountains. Consequently, provinces located in front of these mountain ranges and bordering Laos receive more precipitation than other areas in the region.

5. Conclusion

This study demonstrates that nonstationary GEV models, particularly Model 2, provide a robust framework for forecasting extreme daily precipitation in Thailand. The findings can support sustainable water resource planning and disaster risk management. For future research, the integration of the r-largest method with EVT will be investigated to expand the applicability of extreme value models in real-world hydrological forecasting.