Spatial Modeling of Extreme Temperature in Northeast Thailand

: The objective of the present study was to examine and predict the annual maximum temperature in the northeast of Thailand by using data from 25 stations and employing spatial extreme modeling which is based on max-stable process (MSP) using schlatter’s method. We analyzed extreme temperature data using the MSP using latitude, longitude, and altitude variables. Our result showed that the maximum temperature has an increasing trend. The return level estimates of the study areas from both the local generalized extreme value (GEV) model and MSP models show that the Nong Khai, Maha Sarakham, and Khon Kaen stations had higher return levels than the other stations for every return period, whereas Pak Chong Agromet had the lowest return levels. Furthermore, the results showed that MSP modeling is more suitable than point-wise GEV distribution. We realize that the spatial extreme modeling based on MSP provides more precise and robust return levels as well as some indices of the maximum temperatures for both the observation stations and the locations with no observed data. The results of this study are consistent with those of some previous studies. The increasing trend in return levels could affect agriculture and the surrounding environment in northeast Thailand. Spatial extreme modeling can be beneﬁcial in the impact management and vulnerability assessment under extreme event scenarios caused by climate change.


Introduction
Global climate change has become an urgent matter of concern as it threatens our planet and affects lives and economy both directly and indirectly. Direct effects include, among others, the risk of flooding from rising sea levels, while indirect effects include higher food prices, owing to direct effects on crop cycles and growth. Climate change affects strong weather events, including summer storms, cyclones, and typhoons. It strengthens the El Niño and La Niña events, which are important factors in determining the direction of frequent seasonal storms and the intensification of inclement weather, such as heavy rain, in some countries. The Association of Southeast Asian Nations (ASEAN), including Thailand, is a global food producer that has experienced longer hot and dry spells in their areas of activity, resulting in water and food shortages. According to the Intergovernmental Panel on Climate Change (IPCC) report [1], risks due to climate change are increasing, and should be considered as a part of the 21st century climate system.
Climate change has affected Thailand in several ways; its impacts on the changes in temperature and rainfall are particularly notable. Thus, preparations to deal with and adapt to these changes are essential. In Thailand, when considering the change in temperature since the beginning of record keeping (1951) by the Thai Meteorological Department (TMD) [2], it was found that the average temperature has been increasing. In 2021, the TMD [2] found that the annual average temperature was 27.5°C, which was 0.4°C higher than that in the previous 30-year period (1981-2010). Limsakul [3]'s research supports that day and night temperatures in Thailand have risen considerably and are likely to continue to increase. The northeastern region has shown a more rapid increase in temperature than other areas in Thailand. Increases in temperature impact agriculture, such as through water shortages and limiting plant growth. In addition, the northeastern region frequently experiences severe droughts. Because many people in the northeastern region are farmers, which comprise 46.6% of the population in Thailand [4], rises in temperature, especially during intense heatwaves, will significantly (and likely negatively) influence the means of living of many people in the region.
For the prediction of extreme temperature across Thailand, the statistical models employed by previous studies are rather simple and standard. Sharma and Babel [5] and Manomaiphiboon et al. [6] used regression models and the Mann-Kendall trend detection method. They reported that the extreme temperature tend to increase significantly for all stations examined in Thailand. Limsakul [3] used the empirical orthogonal function approach whereas Rodchuen et al. [7] used an ARMA model. The results of the aforementioned studies indicated that the frequency of hot (cool) temperature extremes has increased (decreased) and will continue in the near future.
Some studies, such as those by [5,[8][9][10], have attempted to employ extreme value models to analyze extreme temperature data in Thailand. Seenoi et al. [10] modeled extreme temperatures in upper northeastern Thailand at nine meteorological stations using a generalized extreme value (GEV) distribution (GEVD). They estimated the parameters using the maximum likelihood estimation (MLE) method under stationary and non-stationary settings. None of the above-mentioned studies applied the model spatially, but only locally; here, local analysis means that the model was built independently for each weather station without considering the spatial dependency among nearby stations. It is generally known that spatial modeling of extreme values can reduce the mean squared error of prediction compared to local modeling [11][12][13]. The reason for the better performance of spatial modeling is that the spatial extreme model takes into account the spatial dependency among nearby stations, whereas the local extreme model is applied independently for each station without considering any events around nearby stations. However, spatial extreme modeling based on the max-stable process (MSP) requires intensive computation; thus, analysis with large number of locations is challenging [11]. Further, the MSP also requires a clear model specification. These disadvantages have prevented researchers from employing the MSP models to analyze extreme spatial data.
In our study, a spatial extreme value model based on an MSP was used and applied to the maximum temperature data in northeastern Thailand. The objective of the spatial analysis is to model a region where extreme events occur, and data are continuously stored [14]. Spatial analysis is an important in-depth analysis, as researchers sometimes cannot determine how independent or close each observation area is. Therefore, the spatial analysis was used to model spatial extreme events. This is a way to analyze spatial data and extreme value properties simultaneously. Yoon et al. [15] studied spatial GEV models within continuous local extreme events using observations of annual daily maximum rainfall data in northeastern Thailand, comprising 25 locations for the period 1982-2013. It was shown that the regional spatial GEV model reflects the spatial pattern well compared to the region-wide spatial GEV model as a local distribution of GEV and provides a strong return level. Several studies have suggested that MSPs are useful for statistical modeling of spatial extremes. These processes are natural extensions of multivariate extreme-value distributions to infinite dimensions [12,[16][17][18]. Thus far, no study has been conducted on the spatial GEVD using observational temperature data in Thailand. Therefore, researcher are interested in studying it to prepare for problems that may arise from temperature changes.
Overall, in this study, we investigated the local GEV and spatial GEVD with extreme values of daily maximum temperature data in northeastern Thailand for 25 stations and predicted the return level of temperature in the 2, 10, 25, 50, and 100 year return periods. Section 2 presents the study areas and meteorology; Section 3 presents the modeling methodology, including the GEV model structure and spatial GEV model; Section 4 describes the results of both local GEV and spatial GEV models; and, finally, Section 5 presents the discussion and summarizes the conclusions of this study.

Study Area
The northeastern region of Thailand, which was considered as the study area, approximately 168,854 km 2 , which is comparable to one-third of the total area of Thailand. This region is 120-400 m above sea level. The northeastern region is mostly a pan-like plateau divided into two large areas: the Korat River Basin, which includes the Mun and Chi rivers, and the Sakon Nakhon River Basin in the northern part of the region.
The general climate in Thailand has an average year-round temperature of 26-27°C, with the lowest average temperature above 18°C; therefore, Thailand is included in the tropical climate zone. The difference in climatic temperatures indicates that the months with the highest air temperature are April to early May, when the temperature reaches 40-41°C, whereas the coldest months are December to January, when the temperature drops below 10°C. The air temperature distribution was found to be high during April-September. When entering October, the air temperature gradually decreases and cools from November to February, and the weather warms again as March begins. From 1951 to 2012, the TMD [2] reported that the annual average temperature for the northeastern region was 28.6°C, but the hottest air temperature was 43.9°C in the Udon Thani Province in late April [19]. Figure 1 shows ridge line plots of the monthly average temperature from 1989 to 2019 for six stations (Nong Khai, Udon Thani, Sakon Nakhon, Khon Kaen, Chaiyaphum, and Nakhon Ratchasima) in the northeast region of Thailand, consistent with the previous description. The northeastern region usually experiences a long period of warm weather because of its inland nature and tropical latitude zone. Between March and May, the peak temperature is about 40°C. In winter, monsoon winds from China cover Thailand, causing occasional cold weather and the temperature to drop to a relatively low value. Particularly in the northeast, the temperatures can drop to near zero°C [19].
In this study, we used yearly data to analyze GEV and then selected the annual maximum temperatures from the daily maximum temperature data obtained from the TMD [2]. The data records were for the duration of 1989 to 2019 from 25 stations in northeastern region of Thailand. The Table 1 shows the station names and identifications, including the altitude above sea level and descriptive statistics of the temperature data (unit:°C) for all 25 stations. All of the stations in the northeastern region of Thailand are shown in Figure 2. Although the weather stations were well distributed, some were located close to each other. Boxplots of the annual maximum temperature data for each station are plotted in Figure 3. This shows that some stations were located close to each other, but had different temperature distributions. For example, the Nakhon Ratchasima (431201), Pak Chong Agromet (431301), and Chok Chai (431401) stations are in the same province but have different temperature characteristics. This may be due to different topographies. The Pak Chong Agromet station is surrounded by mountains, and thus, the ambient temperature is relatively low. The Chok Chai Station is mostly in the highlands and has a river running through it. Therefore, both stations had lower temperatures than the Nakhon Ratchasima station. Table 1. Temperature monitoring station and height above sea level for 25 stations and descriptive statistics of annual maximum temperature data (unit:°C) in northeastern region of Thailand. N and SD represent number of observations and standard deviation, respectively. Q 1 is 25th and Q 3 is 75th percentile.

Methodology
Analysis of the extreme value model with extreme value theory (EVT) can be divided into two types according to the nature of the selection of the extreme value data: GEVD and generalized Pareto distribution (GPD). Modeling with GEVD is suitable for analyzing extreme values over time periods of interest, such as daily, weekly, monthly, quarterly, and yearly. GPD analysis is suitable when there is a large amount of daily data. In this study, we analyzed the GEVD by selecting yearly maximum values from daily temperature data recorded from 1989 to 2019 for 25 stations.

GEVD
The GEVD developed by Jenkinson [20] is flexible to three other distributions: Gumbel, Fréchet, and Weibull distributions. Assuming that X i , where i = 1, 2, · · · , n are independent random variables and have the same probability density function F(x), the maximum value of the random variable is X (n) = Max(X 1 , X 2 , · · · , X n ).

Local GEVD
The cumulative distribution function (CDF) of the GEVD is as follows [11]: where µ, σ, and ξ denote the location, scale, and shape parameters, respectively. The case with ξ → 0 is the Gumbel distribution The cases with ξ > 0 and ξ < 0 are known as Fréchet and negative Weibull distributions, respectively.

MLE
We estimated the parameters using the MLE method. Assuming that X 1 , X 2 , · · · , and, X n are independent variables and have a GEVD, then the log-likelihood function can be written as follows: where ξ = 0 and 1 + ξ x i −µ σ > 0 for i = 1, · · · , n. For the case ξ = 0, it is the limit of the Gumbel distribution; then, the log-likelihood function becomes In practice, it is easier to maximize the log-likelihood function. We obtained the values ofμ,σ, andξ from Equations (3) and (4).
The goodness-of-fit test used for this study was the Kolmogorov-Smirnov (KS) test, which is obtained by transforming to a standard Gumbel distribution, defined by [11].

Return Level Estimation
After the parameters were estimated using the maximum likelihood method, the return level value z p for 0 < p < 1, where z p is defined as the value expected to exceed the average once every 1/p Coles [11], can be calculated as follows: where y p = −log(1 − p). The standard error of the return level was calculated using the delta method. The extreme value model has a preliminary agreement that the distribution of the GEV data at each station is independent of one another; thus, our data are a multivariate series, since data from several locations were recorded. Therefore, the extreme values theory approaches are insufficient and were, thus, modeled using the spatial extreme value, as presented [21].

Spatial GEVD
In this study, we created a spatial model for extreme values via an MSP using Schlater's characterization. The MSP is an infinite-dimensional generalization of multivariate EVT distribution. The concept of spatial dependence is that everything is interrelated, but the nearer ones are more closely related than the distant ones [21], as shown in Figure 4.
For data analysis, we used linear transformations to provide geographic information of the station, as in Equation (7).
to make the inverse calculation of the Hessian matrix more stable [15]. Let M(l, t) be the daily maximum temperature data for the location (l) and period (t) in the spatial domain D ⊂ R 2 [21]. Then, the distribution of M(l, t) is: where µ(l, t), σ(l, t), and ξ(l, t) are the locations, scale, and shape parameters from the GEVD. A stochastic process Z(r) is defined as an MSP if successions of continuous function a n (r) and b n (r) exist, as in Equation (9).
In the present study, we used Schlater's method to rewrite Equation (9) as Equation (10): where U i are the points of a Poisson process at (0, +∞), and Y 1 (x), · · · , Y n (x) are independent replications of the stochastic process Y(x). Portero Serrano et al. [22] presented at bivariate CDF, as follows: where ρ(h) is the correlation function and the distance h separates the different maxima: where τ and η are the scale and shape parameters. The parameter ν represents the nugget effect of measurement errors and microscale variations in the data [13]. The correlation function ρ(h) between two random (X, Y) items should have a positive correlation for ρ = 1 as a full dependency, with this function being constrained by the extremal coefficient (θ) describing the characteristics of the matrix of the dependencies tail and the extremal coefficient function To select the best trend surface model, we used the Takeuchi information criterion (TIC) [25] and, chose the smallest available one. Through a model selection procedure, we built a regression-based form for the location, scale, and shape parameters, as follows: µ(x) = µ 0 + µ lon lon(x) + µ lon lon(x) 2 + µ lat lat(x) + µ lat lat(x) 2 + µ alt alt(x) σ(x) = σ 0 + σ lat lat(x) + σ alt alt(x) (13) ξ(x) = ξ 0 where these models are several combinations of models with the latitude, longitude, and altitude of the station at which the data were observed x. The shape parameter ξ(x) was treated as a constant. The parameters ψ of the model were identified by the pairwise marginal density and estimated by maximizing a composite log-likelihood function, as follows [24]: and TIC function is given by where where i = 1, 2, · · · , n, j = 1, 2, · · · , m and the return level analysis concept are used to estimate the future extreme temperature.

GEVD
From Table 2, the 95% confidence interval estimated the contour shape parameters from the local GEV of the study area, indicating the optimal distribution at each station. There are 10 stations with the contour parameter, and positive coverage means that the tail is skewed to the right; therefore, the optimal distribution is the Gumbel distribution. There are 15 stations in which the shape parameter approximation is in the zero range, that is, the Weibull distribution. We performed our model with the KS test statistic at a significance level of 0.05 by comparing the CDF values of the sample data with the CDF values of the actual data. It was found that the Gumbel and Weibull distributions were appropriate for the data. The return level estimates of the study areas from the local GEVD calculated from Equation (6) are shown in Table 3, indicating that Nong Khai Station (352201) had higher return levels for temperature than the other stations, whereas Surin Station (432201) had the lowest return levels for every return period.

Spatial GEVD
We modeled the spatial model for extreme temperature data for 25 stations in northeast Thailand via an MSP using Schalater's method with a powered exponential covariance function ρ(h) [15,23]. Table 4 shows the parameter estimation with standard error in parenthesis which was calculated via the MSP method referenced in Equation (14). Through a model selection from minimum TIC criterion among possible forms of µ(x) and σ(x) from Equation (15), we obtained the location and scale parameters as in Table 4. From this table, lat 2 represents the squared term of the latitude of x and the shape parameter is treated as a constant. Table 4. The parameter estimation with standard error in parenthesis obtained from spatial GEV models. Calculated via the max-stable process method referenced in Equation (14).  Figure 5 presents the scatter plots of µ, σ, and ξ from the GEVD and MSP for 25 stations, independently corresponding to the model selection procedure. Based on the selected model, we obtained estimates of location (µ), scale (σ), and shape (ξ) parameters for regional spatial GEV in Table 5.   Figure 6 show the return levels corresponding to 2, 25, 50, and 100 year return periods, which were obtained from an MSP using the regional models for 25 stations in the northeast region of Thailand. Table 6 shows that the stations with the highest return levels were Nong Khai (352201), Maha Sarakham (387401) and Khon Kaen (381201), while Pak Chong Agromet station (431301) had the lowest return temperature levels. From Figure 5. Relationships between geographic covariance and GEV parameters where parameter estimates were obtained from GEV models for 25 stations independently. Table 5. Estimated parameters and standard error values (in parenthesis) obtained from regional spatial GEV model.

Station
Locationμ ( Table 6 and Figure 6 show the return levels corresponding to 2, 25, 50, and 100 year return periods, which were obtained from an MSP using the regional models for 25 stations in the northeast region of Thailand. Table 6 shows that the stations with the highest return levels were Nong Khai (352201), Maha Sarakham (387401) and Khon Kaen (381201), while Pak Chong Agromet station (431301) had the lowest return temperature levels. From Figure 6, we can see that the MSP model collects geographical and covariate information well across the region. We used a Kriging and inverse distance weighting (IDW) technique for interpolation in drawing Figure 6. Then we compared the IDW and Kriging techniques and found that the IDW method was more effective. See detailed results in supplementary material. This result is consistent with the TMD report that these stations are relatively flat and the climate is classified as tropical [2].

Discussion and Conclusions
In this study, we used GEV and spatial GEV modeling using the MSP to analyze the annual maximum temperatures in northeastern Thailand at 25 stations. We used Schlathers method with a powered exponential correlation function to build a regional model. In addition, we predicted the return levels of the observation data for several return periods using regional models and drew return level maps. This study showed that 15 stations had Weibull distributions and 10 had Gumbel distributions, which were suitable for the significant KS tests for all stations for the local GEV model. The return level estimates of the study areas from the local GEV showed that Nong Khai (352201) had higher return levels than the other stations for every return period. For the spatial GEV model, we analyzed extreme data using MSP with latitude (lat), longitude (lon), and height above sea level (altitude, alt) variables. The equations suitable for modeling are as follows: We set the shape parameter as constant for the spatial GEV models. Thus, future research could consider a nonconstant shape parameter ξ(x), which may require more samples. Building the best regression-based form for the location and scale parameters requires a substantial computational burden. We selected covariates at the minimum TIC from many possible models. We employed a stage-wise selection procedure to select the best regressionbased form from all the possible combinations. However, we observed failure of the TIC computation in a few cases, which made the procedure unstable. In this regard, future research should develop a fast and stable algorithm for the better selection of variables. The Figure 6 shows the MSP model of the spatial GEV, which collects geographical and covariate information across the region. The stations with the highest return temperatures were Nong Khai (352201), Maha Sarakham (387401), and Khon Kaen (381201), whereas Pak Chong Agromet (431301) had the lowest return levels in every return period. The estimation results of the return level of our study is consistent with Limsakul's [3] research showing that Thailand's temperature tends to increase, as well as Seenoi's [10] research supporting that temperatures in upper northeastern Thailand tend to increase.
Due to the continually increasing estimate of the return level, the northeastern region of the country may be affected by some risky extreme events. These are, for example, severe droughts and intense heatwaves that could result in forest fires or significantly fewer agricultural products. This is coupled with the possibility of overall temperatures exceeding 35°C and the spread of pests and plant pathogens.
The approach of this study can be applied to the field of numerical model outputs of climate systems and compared to a model fitted to climate observations. Future studies should seek a significant covariate that affects temperature. The latitude, longitude, and altitude data are considered important; however, we may have missed more important covariates for extreme temperatures in this study, such as topographic aspects and coastal proximity. One of the reasons why the global spatial GEV model is not a good reflection is that the geographic covariates do not reflect the characteristics of the region. If the covariates that have influence are considered in the model, it is expected that the explanatory power of the model will be further enhanced. Furthermore, the detection and prediction of changes in climate, including trends in extreme temperatures [13] and extreme wind speeds [26], will be performed. In addition, the MSP modeling approach for different time periods can be used to determine climate change extremes in Thailand.
Modeling the entire distribution (including the maximum, mean, and minimum simultaneously) of a climate variable for a spatial field and using it to detect changes (e.g., shifts in means and standard deviations) may be more reliable and valuable than modeling only extremes [27]. For this purpose, a time-dependent MSP model with climatic covariates [28] as well as a spatial-temporal linear model [26,29] can be useful. Addressing these challenges may be the first step in delving further into this research area. Moreover, increased collaboration between climate scientists and statisticians is required.