On the use of an innovative trend analysis methodology for temporal trend identification in extreme rainfall indices over the Central Highlands, Vietnam

Understanding past changes in the characteristics of climate extremes forms an essential part of viable countermeasures to cope with climate-induced risks in a rapidly warming world. Thus, this paper endeavored to explore possible non-monotonic trend components in heavy rainfall events over the Central Highlands (Vietnam) by employing Şen’s innovative trend analysis method in conjunction with the well-defined extreme rainfall indices. The outcomes show less spatially coherent trends in the intensity, frequency, and duration of extreme rainfall events across the study area, and most analyzed indices exhibited non-monotonic trend forms. The overall trends in the intensity and frequency, represented by maximum 1-day, 5-day precipitation amount (Rx1day, Rx5day), very and extremely wet days (R95p and R99p), and the number of very and extremely heavy precipitation days (R20mm and R50mm), were characterized by significant increases at three or four sites. Concerning dry and wet extremes, observed increases in consecutive dry days (CDD) and decreases in consecutive wet days (CWD) were identified significantly simultaneously at three sites. There is high domination of significant increases in high-value subgroups of these indices. Several sites also exposed significantly increasing trend behaviors in these indices within all low-, medium-, and high-value subgroups, thereby implying the intensity and frequency have intensified over recent decades. This study also indicated potential linkages between annual peaks of several extreme rainfall indices, characterizing intensity, frequency, and wet duration, usually coincided with negative Indian Ocean Dipole and La Niña events, while dry extremes were usually found during El Niño events based on correlation analysis.


Introduction
The vast majority of anthropogenic greenhouse gas emissions from fossil fuel combustion have been releasing into the Earth's atmosphere since the mid-nineteenth century, which gives rise to profound changes in the Earth's climate.
An up-to-the-minute assessment report released by the Intergovernmental Panel on Climate Change (IPCC) documented that the global surface temperature for the most recent decade (2011-2020) was 1.09 [0.95 to 1.20] °C higher than the average over 1850-1900 period, and also projected to rise by 1.2-1.9 °C in the near-term 2021-2040, by 1.2-3.0 °C in the mid-term 2041-2060, and by 1.0-5.7 °C in the longterm 2081-2100 compared to the average of the period 1850-1900 under five new illustrative emissions scenarios (IPCC 2021). Moreover, human beings experience both changes in mean climate and extreme events such as hot/cold days/nights and heavy precipitation events. These moderate and rare weather and climate extremes pose a significant threat to society and biophysical systems (Dunn et al. 2020;World Meteorological Organization 2009;Zhang and Zwiers 2013). Thus, it is imperative to address trend possibilities of weather and climate extremes in the context of a warmer climate. It is advisable to utilize appropriate descriptive indices of extremes and extreme-value theory (World Meteorological Organization 2009).
From a global perspective on analyzing weather and climate extremes, a core set of 27 extremes indices for temperature and precipitation (World Meteorological Organization 2009;Zhang et al. 2011;Zhang and Zwiers 2013), defined by the Joint CCl/CLIVAR/JCOMM Expert Team on Climate Change Detection and Indices (ETCCDI) is adopted worldwide. Zhang et al. (2011) documented various evident indications of increasing (decreasing) trends in warm (cold) extremes across the globe, while heavy precipitation events tend to be more frequent and intense, which is attributable to human-induced increases in greenhouse gas concentrations (Min et al. 2011). Donat et al. (2013) and Dunn et al. (2020) developed the well-known global gridded datasets of climate extremes, namely GHCNDEX and HadEX3, respectively, which yield numerous applications in extreme analysis (such as trend identification, model evaluation, climate change detection, and attribution). Dunn et al. (2020) also substantiated several discrepancies when using 1961-1990 versus 1981-2010 reference periods for calculating threshold-based indices.
One of the most common themes is comparative investigations into the performance of diverse data sources in current and future assessments of extreme events. Alexander et al. (2020) conducted an intercomparison between in situ, satellite, and reanalysis products in reproducing global precipitation extremes and highlighted the applicability of satellite-derived precipitation datasets, especially in ungauged or poorly gauged regions. Kim et al. (2020) demonstrated limited improvements in the Coupled Model Intercomparison Project phase 6 (CMIP6) model skills for the 27 ETCCDI climate extremes indices compared to the CMIP5 models. Lately, Sun et al. (2021) indicated significant increases in the annual maximum one-day (Rx1day) and five-day (Rx5day) precipitation at the global and continental scales, as well as in many large regions over the periods 1900-2018 and 1950-2018, thereby implying the observed intensification of extreme precipitation. Based on the state-of-the-art CMIP6 simulations under four socioeconomic development pathways, Li et al. (2021) substantiated that the frequency and intensity of hot temperature and precipitation extremes are likely to increase significantly over most of the Earth's surface.
Another crucial topic under consideration is the identification of long-term temporal trends in the ETCCDI climate extremes indices. Historical changes in the extreme aspects of climate (such as intensity, duration, and frequency) play an essential role in projecting future climate change scenarios (Min et al. 2011), as well as counteracting climaterelated risks to humans, ecosystems, and infrastructure (World Meteorological Organization 2009). Over the last few decades, there has been a proliferation of case studies that employ parametric and non-parametric statistical trend tests for analyzing past changes in extremes indices. Dimri (2019) applied the modified Mann-Kendall (MK) test and Sen's slope estimator to compare regional and seasonal trends in eight extreme temperature indices derived from the ERA-Interim and IDM datasets over India. Dookie et al. (2019) utilized the least-squares regression and Student's t-test to demonstrate local patterns in eleven temperature and ten precipitation indices over two Caribbean small islands. García-Cueto et al. (2019) applied the MK and Sen's slope tests in 14 ETCCDI climate extremes indices to indicate changing climate conditions over some urban areas in Mexico. Murara et al. (2019) also employed the linear least square method and MK test to estimate the magnitude and statistical significance of trends in ten extreme precipitation indices in the Itajaí river basin (Brazil). Qian et al. (2019) applied various parametric and non-parametric methods (such as ordinary and refined least squares regression, classical and modified MK, and Sen's slope) in eight extreme temperature indices derived from observational and gridded datasets across China. Qian et al. (2019) also emphasized the prime importance of non-Gaussian and serial independent assumptions when applying classical parametric and non-parametric methods to analyze trends in climate extremes indices. Costa et al. (2020) employed linear regression and sequential MK test to estimate annual trends in ETCCDI climate extremes indices over northeast Brazil. Venkata Rao et al. (2020) explored spatio-temporal changes in nine extreme rainfall indices by utilizing Sen's slope, classical MK, and its modifications over two floodprone basins in eastern India at three time periods (i.e., longterm (1901-2018), pre-1950, and post-1950).
In Vietnam, the characteristics of climate extremes are assessed comprehensively at the national and subregional scales. Ho et al. (2011) delved into the past and future trends in hot days, cold nights, and heavy rainfall days over seven climatic subregions in Vietnam based on observed records and Regional Climate Model (RCM) outputs. Ngo-Duc et al. (2014) further applied Student's paired t-test, Sen's slope, and MK test to explore projected changes in the annual maximum value of daily maximum temperature (TXx), the annual minimum value of daily minimum temperature (TNn), and Rx1day over Vietnam based on the outputs of three RCMs. At the city scale, Quan et al. (2021) presented a well-conducted investigation into the statistical strength, stability, and magnitude of historical trends in seven extreme rainfall indices in Ho Chi Minh City (Vietnam) by employing Sen's slope and MK test in conjunction with the trendfree pre-whitening procedure. Khoi et al. (2021) further applied the same trend analysis methods to explore future changes in the spatio-temporal patterns of extreme rainfall indices over Ho Chi Minh City (Vietnam) based on the statistically downscaled climate data.
It is acknowledged that most previous trend search studies on climate extremes were carried out mainly through parametric and non-parametric approaches such as linear regression, Sen's slope estimator, classical and modified MK tests. These commonly used statistical methods need to comply with some assumptions such as normal distribution, serial independence, and sufficient length of observed records. In practice, almost all hydro-meteoro-climatological data are serially correlated, which will lead to misinterpretation of trend existence in a given time series (Almazroui and Şen 2020;von Storch 1995). Qian et al. (2019) also indicated that most climate extremes indices are unlikely to follow a Gaussian distribution as well as serial independence. Additionally, these aforementioned methods identify only overall linear or monotonic trends within a given considered period. On the other hand, an innovative trend analysis (ITA) methodology, originated by Şen (2012, 2014, 2017), can be employed to identify temporal trends in a given time series without considering any restrictive assumptions (Almazroui and Şen 2020; Şen 2012). This cutting-edge method also yields more detailed information about trend behaviors such as non-monotonic trends and subtrends in a given time series (Dabanlı et al. 2016;Elouissi et al. 2016;Şen 2012). Recently, the ITA method was applied successfully to identify trend possibilities of conventional hydro-meteorological variables (Dabanlı et al. 2016;Elouissi et al. 2016;Öztopal and Şen 2017), as well as derived indices such as Standardized Precipitation Index (Caloiero 2018;Danandeh Mehr and Vaheddoost 2020;Phuong et al. 2021;Yilmaz 2019), Standardized Precipitation Evapotranspiration Index (Danandeh Mehr and Vaheddoost 2020), and Effective Drought Index (Malik et al. 2020).
In general, applying the ITA method in climate extremes indices can yield more insights into the changing character of climate extremes. Thus, the present study aimed to delve deeper into the overall and partial trends in extreme rainfall events over the Central Highlands of Vietnam by employing the ITA method in conjunction with ten ETCCDI extreme rainfall indices. This study also addressed potential linkages between interannual variability of extreme rainfall events and large-scale drivers such as El Niño-Southern Oscillation (ENSO) and Indian Ocean Dipole (IOD) based on correlation analysis. It is expected that understanding temporal trends in low-and high-value subgroups of extreme rainfall indices could contribute more scientific merits to local decision-makers and other stakeholders to cope with climaterelated risks under a changing climate.

Study area
The selected study area is the Central Highlands of Vietnam ( Fig. 1), which is a mountainous region with elevation ranging considerably between 100 and 2400 m a.s.l.. The Central Highlands, covering an area of approximately 54,508 km 2 (accounting for around 16.5% of the total area of Vietnam), lies approximately between the latitudes 11° 20′ N-15° 40′ N, and the longitudes 107° 20′ E-109° 00′ E. In 2020, the number of inhabitants was around 5.9 million persons (accounting for nearly 6.1% of the total population of Vietnam). Diverse ethnic groups are living in the Central Highlands, and the socioeconomic conditions are relatively low compared to other regions in Vietnam.
There are 12 meteorological stations located quite proportionally to the whole extent of the Central Highlands ( Fig. 1). Table 1 shows some basic information about geographical location and several statistical characteristics of annual total rainfall at these sites. It is noting that the available rainfall records cover 38-year (1982-2019) to 59-year (1961-2019) periods depending on certain meteorological stations. According to the well-known world map of the Köppen-Geiger climate classification proposed by Peel et al. (2007), the general climate of the Central Highlands is categorized mainly as tropical savannah (Aw) climate. The average annual temperatures in the Central Highlands are approximately 18.0 °C-25.9 °C, which is significantly lower than the temperatures in surrounding climatic subregions (Ho et al. 2011). Additionally, there are two distinct seasons, namely dry and rainy seasons, in the Central Highlands. Rainfall during the rainy season (May-October) contributes mainly to the total amount of annual rainfall, varying from around 1200 to 2700 mm. Ho et al. (2011) documented a coincidence of the rainy season in the Central Highlands and the active phase of the Asian monsoon.

Extreme precipitation indices
Moderate extremes typically occur at least once or several times every year (World Meteorological Organization 2009; Zhang et al. 2011). It is recommended to employ an internationally coordinated core set of 27 ETCCDI indices to comprehend potential changes in the characteristics of moderate   Zhang et al. 2011;Zhang and Zwiers 2013). The ETCCDI indices, derived from daily weather data, were defined to facilitate many difficulties in understanding past, current, and future changes in weather and climate extremes. Additionally, utilizing agreed-upon indices makes it possible to investigate comparative studies of the characteristics of extremes in different climatic regions across the globe and draw out a global perspective on climate extremes accordingly .
In essence, the ETCCDI indices can be classified into three types Zhang and Zwiers 2013). The first one is absolute extremes (e.g., TXx, TXn, TNx, TNn, Rx1day, or Rx5day), quantifying the monthly/annual maximum/minimum values of daily temperature and maximum values of daily precipitation amounts, which is substantially conducive to a wide range of engineering fields. The second and third types are day-count indices with fixed and percentile thresholds, respectively. The day-count indices with fixed thresholds (e.g., SU, TR, R10mm, R20mm, or R50mm) are more suitable for climatic-change impact studies on health, hydrology, and agriculture. Meanwhile, the day-count indices based on percentile thresholds (e.g., TX10p, TX90p, TN10p, TN90p, R95p, or R99p), expressing anomalies relative to the local climate, allow for meaningful spatial comparisons of different regions with different climate conditions. The present study analyzed ten extreme precipitation indices as outlined in Table 2. These indices were calculated on an annual basis. It is noting that R50mm is a user-defined index. The user-defined threshold was set to 50 mm since this extreme threshold is currently used by the Vietnam Meteorological and Hydrological Administration. The R50mm was also analyzed by Endo et al. (2009); Ngo-Duc et al. (2014); Tangang et al. (2018). For the percentile-based indices, it is necessary to determine a reference period explicitly. Based on the availability of meteorological datasets, this study calculated percentile-based indices using a 1981-2010 reference period, which is recommended by World Meteorological Organization (2017) and used by Dunn et al. (2020). This study utilized the R programming language (Ihaka and Gentleman 1996;R Core Team 2021) in conjunction with the "climdex.pcic" package (David Bronaugh for the Pacific Climate Impacts Consortium 2020) for extreme indices calculation and the "tidyverse" package (Wickham et al. 2019) for data manipulation and visualization. Şen (2012, 2014, 2017) developed an innovative trend analysis methodology, which is fundamentally based on constructing a scatter plot of a set of subseries derived from a given time series on a Cartesian coordinate system. Furthermore, various improved versions of the ITA method were proposed to enhance visual properties and assess statistical significance (Phuong et al. 2021). In comparison with the well-known (non)parametric statistical trend tests such as linear regression and MK methods, the ITA method does not require any restrictive assumptions such as normality and serial independence (Şen 2012, 2014). Another advantage of the ITA method is that this method can be applied to a small sample (e.g., ten observations), while the other classical approaches necessitate at least 30 data records (Almazroui and Şen 2020). Especially, non-monotonic trends (i.e., the composition of various trend components) in a given hydro-meteorological time series can be detected by the ITA method (Şen 2012, 2014). Specifically, the ITA method can yield more detailed information about the overall and partial trends based on deliberate categorization. Figure 2 depicts a typical innovative trend plot. A step-by-step procedure to implement the ITA method  (2018) and Öztopal and Şen (2017). In particular, the first step is to divide the original time series into two (or more) equal subseries and then sort each subseries in ascending order after preparing well-founded datasets. Then, the first half subseries is plotted on the abscissa axis against the second one on the ordinate axis. The final step is to draw the 1:1 (45°) straight line and ± 10% lines on the Cartesian coordinate system. Theoretically, the 1:1 (45°) straight line divides the ITA plot into two equal regions (Fig. 2). The upper/ lower triangular areas indicate increasing/decreasing trends, while the 1:1 (45°) straight line corresponds to trend-free cases. Additionally, Şen (2012, 2014) substantiated that the trend slope tends to get stronger when the scatter points fall far from the 1:1 (45°) straight line. To yield more quantitative interpretations, one can calculate the trend slope developed by Şen (2017) according to the following expression:

Innovative trend analysis methodology
where s is the ITA trend slope. n is the sample length. y 1 and y 2 are the arithmetic averages of the first and second halves of a given time series, respectively. s = 2 y 2 − y 1 n 3 Results

Overall and partial trends in extreme rainfall indices
Temporal trend possibilities of ten ETCCDI extreme rainfall indices over the Central Highlands of Vietnam were explored by employing the ITA method, as illustrated in Figs. 3,4,5,6,7,8,9,10,11,and 12. In these figures, the blue diamond symbols are centroid points used for indicating overall trends, while the gray-filled circle, square, and triangle symbols are mean points of low-, medium-, and high-value clusters, respectively, used for identifying corresponding subtrends as suggested by Dabanlı et al. (2016). Figures 3 and 4 depict the results of overall and partial trend identification of the annual maximum amount of precipitation accumulated over one day (Rx1day) and five consecutive days (Rx5day), respectively. The overall trends in Rx1day and Rx5day have neutral behaviors at most stations, accounting for 7 out of 12 stations. The increasing trends in Rx1day were obtained in four stations (i.e., Kon Tum, Buon Me Thuot, Dak Nong, and Bao Loc), with the estimated ITA trend slope varying around 0.73-0.89 mm/year, while the Ayun Pa station was the only one responsible for the significant declining trend at the rate of approximately − 1.74 mm/year. In comparison, three stations (i.e., Kon Tum, Buon Me Thuot, and Bao Loc) also experienced significant increases in Rx5day at the rate of around 1.17-2.38 mm/year, while significant decreases were found in Dak To and Ayun Pa stations, with the estimated ITA trend slope of approximately − 1.5 mm/ year. Moreover, Figs. 3 and 4 delineate several partial trend patterns in Rx1day and Rx5day. The low-and mediumvalue subgroups were characterized mainly by trend-free behaviors at most stations, accounting for 8-10 out of 12 stations. Concerning the high-value cluster, four stations (i.e., An Khe, Buon Me Thuot, Dak Nong, and Bao Loc) exhibited significant increases in Rx1day, and six stations (i.e., Kon Tum, Buon Ho, Buon Me Thuot, Dak Nong, Da Lat, and Bao Loc) exposed significant increases in Rx5day. It is noting that the Kon Tum station experienced significant increases in Rx5day within all subgroups. The findings of significant increases in the high-value subgroups indicated that the absolute extremes have intensified over recent decades whether the overall trends or the low-and medium-value clusters were significant or not.
Turning to the simple daily intensity index (SDII), Fig. 5 shows that more than half of the analyzed stations exhibited significant increases (including Kon Tum, Pleiku, Buon Ho, M'Drak, Buon Me Thuot, Dak Nong, and Bao Loc). The magnitude of these upward trends was estimated approximately from 0.04 to 0.17 mm/day/year. Concerning subtrend interpretations, it is discernible that the majority of scatter points in high-value SDII subgroups at four stations (i.e., Buon Ho, M'Drak, Dak Nong, and Bao Loc) fall father to the 1:1 (45°) straight line compared to the low-and medium-value clusters, thereby indicating higher increasing trends in high-value subgroups. Moreover, the Bao Loc station exhibited increasing trend behaviors within all subgroups.
Figures 6 and 7 represent temporal trend possibilities in the very wet days (R95p) and extremely wet days (R99p), respectively. As expected, the R95p and R99p exposed the same overall and partial trend patterns over the study area. In particular, 7 out of 12 stations exhibited almost trendless behaviors. Meanwhile, significant increases were found in three stations (i.e., Kon Tum, Buon Me Thuot, and Dak Nong), with the estimated ITA trend slopes varying from 4.3 to 9.4 mm/year as for R95p and from 2.3 to 6.1 mm/year as for R99p. Conversely, significant decreases at the rate of approximately from − 5.0 to − 6.5 mm/year as for R95p and − 2.8 to − 4.7 mm/year as for R99p were obtained in two stations (i.e., Dak To and Ayun Pa).
Figures 6 and 7 further portray potential trends in low-, medium-, and high-value R95p and R99p subgroups. In particular, there is high domination of neutral trends in low-and medium-value clusters as for both R95p and R99p. Meanwhile, the high-value subgroups were dominated mainly by significant upward trends. Eight stations (i.e., Kon Tum, An Khe, Buon Ho, M'Drak, Buon Me Thuot, Dak Nong, Da Lat, and Bao Loc) experienced significant increases in R95p, and six stations (i.e., Kon Tum, An Khe, M'Drak, Buon Me Thuot, Dak Nong, and Bao Loc) exhibited significant increases in R99p. Moreover, the Kon Tum station was the only one responsible for significant increases in R95p and R99p within all clusters. In general, the vast majority of significant increasing trends in high-value subgroups imply that these percentile-based indices have become more intense over recent decades. Turning to the number of days per year with rainfall amount ≥ 10 mm (R10mm), ≥ 20 mm (R20mm), and ≥ 50 mm (R50mm), Figs. 8, 9, and 10 show that few stations exhibited significant increases. In particular, the overall increasing trends in R20mm were found in three stations (i.e., Kon Tum, An Khe, and M'Drak), with the estimated ITA slopes ranging approximately from 0.14 to 0.24 days/ year. Similarly, three stations (i.e., Kon Tum, Buon Me Thuot, and Dak Nong) exhibited significant increases in R50mm at the rate of around 0.05-0.10 days/year. Meanwhile, there is high domination of neutral and decreasing trends in R10mm. Figures 8, 9, and 10 also reveal that the low-and medium-value subgroups of R10mm, R20mm, and R50mm were characterized mainly by trendless and decreasing behaviors. Concerning the high-value subgroups, there were two stations (i.e., An Khe and M'Drak) exhibiting significant increases in R10mm, and five stations (i.e., Kon Tum, An Khe, M'Drak, Dak Nong, and Bao Loc) experiencing significant increases in R20mm. Moreover, most stations exposed significant upward trends in R50mm, accounting for 8 out of 12 stations such as Kon Tum, An Khe, Buon Ho, M'Drak, Buon Me Thuot, Dak Nong, Da Lat, and Bao Loc. It is worth noting that the Kon Tum station demonstrated increasing trend behaviors irrespective of low-, medium-, or high-value subgroups. In general, these outcomes indicate that the very heavy rainfall days (R20mm) and extremely heavy rainfall days (R50mm) have occurred more frequently during the study period. Figures 11 and 12 describe possible trend components in the consecutive dry days (CDD) and the consecutive wet days (CWD), respectively. In particular, significant upward trends in CDD were obtained in three stations (i.e., Ayun Pa, Buon Ho, and Bao Loc) with the estimated ITA trend slopes being around 0.66-1.27 days/year, while two stations (i.e., Kon Tum and An Khe) experienced significant declining trends at the rate of approximately from − 0.23 to − 0.7 days/year. Meanwhile, the overall trends in CWD were characterized by significant declines at four stations Concerning partial trend interpretations of CDD and CWD, Figs. 11 and 12 further show that the low-and medium-value clusters were dominated mainly by trend-free behaviors. Meanwhile, six stations (i.e., Dak To, Ayun Pa, Buon Ho, M'Drak, Lien Khuong, and Bao Loc) exposed significant increasing trends in the high-value subgroups of CDD. Two stations (i.e., Ayun Pa and Bao Loc) also exhibited significant decreasing trends in the high-value subgroups of CWD. It is worth noting that the Ayun Pa station was the only one responsible for significant increasing trends in CDD, and significant decreasing trends in CWD within all low-, medium-, and high-value subgroups. On the whole, these findings indicate that dry conditions in the study area have become more frequent and lasted longer during the analyzed period.

Rx1day | Dak To
It is worth pointing out that there is much heterogeneity of significant changes in extreme rainfall indices among 12 analyzed sites located in the Central Highlands of Vietnam.
For instance, the overall and partial trends at the Kon Tum station, located in the northern part of the Central Highlands, were characterized mainly by significant increases in the extreme aspects of intensity and frequency, represented by Rx1day, Rx5day, SDII, R95p, R99p, R20mm, and R50mm. Moreover, the Kon Tum station was the only one responsible for significant increases in all low-, medium-, and high-value subgroups of Rx5day, R95p, R99p, R20mm, and R50mm. Meanwhile, the Buon Ma Thuot, Dak Nong, and Bao Loc stations are those located in the central and southern parts of the Central Highlands, respectively, exhibiting significant increases in Rx1day, Rx5day, SDII, R95p, R99p, and R50mm; but unlike Kon Tum station, none of these sites exposed significant increases in these extreme indices within all clusters.
By contrast, the Dak To and Ayun Pa stations, also situated in the northern part of the Central Highlands, were dominated mainly by significant declining trends in the extreme aspects of intensity and frequency, represented by Rx1day, Rx5day, R95p, R99p, R20mm, and R50mm. Concerning the dry and wet duration, the Ayun Pa, Buon Ho, Dak Nong, and Bao Loc stations are those located in northern, central, and southern parts of the Central Highlands, respectively, exhibiting significant increases (decreases) in CDD (CWD) simultaneously. Moreover, the Ayun Pa station was the only one responsible for significant declines in R20mm, R50mm, and CWD and significant increases in CDD irrespective of low-, medium-, or high-value subgroups. The total amount of annual rainfall at the Ayun Pa station is substantially lower than the figures at the remaining sites (see Fig. 1 and Table 1) and also experienced a significant reduction over recent decades, as documented by Phuong et al. (2021). It is noting that the Ayun Pa station is located in the northern part of the study area, with the lowest elevation compared to the other sites.

SDII | Dak To
Additionally, it is not uncommon for nearby sites to exhibit opposite trends in the same precipitation-related indices. The most striking case is the contrastive outcomes of overall and partial trends in Rx5day, R95p, R99p, and R50mm at Kon Tum and Dak To stations, which are established in the same province located farther north of the Central Highlands, thereby suggesting less spatial coherence on the whole. Caesar et al. (2011), Choi et al. (2009), and Powell and Keim (2015 also emphasized the low spatial coherence in temporal trend possibilities of ETCCDI extreme precipitation indices.

Links to large-scale drivers
This study further investigated the potential linkages between extreme rainfall events over the Central Highlands of Vietnam and large-scale ocean-atmospheric circulation such as El Niño-Southern Oscillation (ENSO) and Indian Ocean Dipole (IOD). Figure 13 presents the estimated values of Spearman's rank correlation coefficient, which assess the relationships between ten ETCCDI extreme rainfall indices at 12 considered sites, and the Niño 3.4 and dipole mode index (DMI), represented the impact of anomalies in sea surface temperature (SST) in the Pacific and the Indian Ocean, respectively. The Niño 3.4 index, defined as the average of SST anomalies over the region 5° N-5° S and 170°-120° W, is greatly encouraged to monitor and characterize ENSO phases in the Pacific Ocean (Bamston et al. 1997). Meanwhile, the DMI describing the difference in SST anomalies between the western equatorial Indian Ocean (50° E-70° E and 10° S-10° N) and the southeastern equatorial Indian Ocean (90° E-110° E and 10° S-0° N), can be utilized to represent the dipole mode in SST in the Indian Ocean region (Saji et al. 1999). It is noting the Niño 3.4 and the DMI datasets are provided by the NOAA Physical Sciences Laboratory, Boulder, Colorado, USA, from their website at https:// psl. noaa. gov/ gcos_ wgsp/ Times eries/. Figure 13 reveals weak and moderate correlations of the Niño 3.4 index and DMI with the extreme rainfall indices at most sites in the study area, whether the correlations are positive or negative. The outputs with significant correlations with the Niño 3.4 index and DMI are 20 and 14, accounting for approximately 16.67% and 11.67% of the total number of correlation calculations, respectively. Concerning potential linkages between large-scale drivers and the intensity of rainfall extremes, the Rx1day and Rx5day exhibited significant negative correlations with the DMI at An Khe, Ayun Pa, and Buon Ho stations, and the annual highest and lowest values of Rx1day and Rx5day at these sites were usually coincident with negative IOD (e.g., 1992 and 2016) and positive IOD (e.g., 1982IOD (e.g., , 1983IOD (e.g., , 1994IOD (e.g., , 2015IOD (e.g., , and 2018, respectively. The outputs also indicate significant positive and negative correlations between the SDII and the Niño 3.4 index at the Buon Ho and Kon Tum stations, respectively, while a significant positive correlation of the SDII with the DMI was found only at the Bao Loc station. Meanwhile, there were few significant negative correlations of R95p and R99p with the Niño 3.4 index (at Kon Tum and M'Drak stations) and DMI (at Ayun Pa, Buon Ho, and Buon Me Thuot stations), and the annual peaks of R95p and R99p at these sites were usually coincident with negative IOD (e.g., 1992) and strong La Niña event (e.g., 2008).
The correlation analysis also indicated that the frequency of rainfall extremes, represented by R10mm, R20mm, and R50mm, was negatively correlated with the Niño 3.4 index and DMI at most sites. In particular, the annual peaks of R10mm, R20mm, and R50mm were found mainly during negative IOD (e.g., 1974, 1981and 1996) and La Niña events (e.g., 1996, 1998-2001, 2006-2008, 2010-2017. Concerning interannual variability of dry and wet extremes, represented by the CDD and CWD, in relationship with large-scale drivers, positive correlations of CDD with the Niño 3.4 index were found to be statistically significant at four sites (i.e., Ayun Pa, Buon Ho, M'Drak, and Buon Me Thuot). The annual values of CDD at these sites varied approximately from 42 to 100 days during analyzed periods, with the peaks reaching mainly from 127 to 158 days. Moreover, the highest values of CDD at these sites usually coincided with El Niño events (e.g., 1983, 1994, 2003-2005, 2014-2016, 2018-2019), including some very strong El Niño events such as during 1983 and 2015-2016. Meanwhile, the CWD showed a statistically negative correlation with the Niño 3.4 index at the M'Drak station.
The annual values of CWD at the M'Drak station altered mainly from 10-16 days, and its peaks reaching around from 22-30 days, were found to be coincident with moderate and strong La Niña events (e.g., 1995-1996, 1999-2000, and 2011). In addition, the CWD at the Ayun Pa station, varying mainly from 7 to 14 days, was negatively correlated with the DMI. The annual peaks of CWD at the Ayun Pa station also coincided with negative IOD (e.g., 1998).

Discussion
Taking all aspects into consideration, it is evident to conclude that the frequency and intensity of extreme rainfall events in the Central Highlands of Vietnam have increased substantially, which is reflected by the domination of significant increases in high-value subgroups of most analyzed indices (such as Rx5day, SDII, R95p, R99p, R50mm, and CDD). These findings are highly consistent with previous in-depth studies over Southeast Asia and the Indo-Pacific region (Caesar et al. 2011;Endo et al. 2009). In particular, Endo et al. (2009)  heavy precipitation increased (decreased) over the southern (northern) parts of Vietnam during the period 1950s to 2000s. Endo et al. (2009) also revealed significant increases in SDII, R50mm, and CDD over the Central Highlands of Vietnam, which is highly analogous to the above-mentioned findings. Caesar et al. (2011) demonstrated an increasing trend in R95p, while there was no significant change in Rx5day over the Indo-Pacific region from 1971 to 2005. Caesar et al. (2011) also pointed out several increases (decreases) in Rx5day over the southern (northern) parts of Vietnam, which is partly consistent with those in the Central Highlands as described above. Furthermore, it is expected that the extreme rainfall events in the Central Highlands of Vietnam are likely to become more frequent and intense in the future. Tangang et al. (2018) performed multi-model simulations and projections of four extreme precipitation indices over Southeast Asia under global warming of 2 °C and indicated significant increases in R50mm and Rx1day in the Central Highlands of Vietnam. Phuong et al. (2021) pointed out a potential reduction in the total amount of annual rainfall and an increase in the highest intensity of meteorological drought in the Central Highlands by analyzing the lowest values of the 3-month and 12-month Standardized Precipitation Index in each year for the considered period. The present study further substantiated that half of the analyzed stations exhibited significant increases in the high-value subgroups of CDD. These outcomes indicate that the Central Highlands of Vietnam experienced a longer duration of dry conditions during recent decades, especially in the dry season, which could lead to a wide range of constraints on agricultural and other socioeconomic activities in the study area.
This study evaluated moderate extremes by analyzing the well-defined ETCCDI extreme precipitation indices. Further investigations should delve deeper into the intensity and frequency of rare extreme events, which might occur once in 20, 50, or 100 years. To do this, one can employ extreme-value theory such as the "peaks-over-threshold" or "block maximum" method as recommended by World Meteorological Organization (2009). These approaches are usually applied to the absolute extremes (i.e., Rx1day and Rx5day). Li et al. (2021) highlighted the prime importance of understanding changes in the frequency and intensity of rare extremes because such events often cause more adverse impacts on the society and biophysical systems, and most infrastructure design requires estimation of rare extreme events as well. Moreover, human-induced global warming is bound to cause stronger influences on rare extreme events than moderate ones. Apart from analyzing individual indices, it is advisable to examine changes in the ratio of a given pair of indices (e.g., the percentage of annual precipitation due to very or extremely wet days). Such an analysis yields more understanding of relatively larger or smaller changes in extreme precipitation events than the total amount (World Meteorological Organization 2009). Another important theme is to explore potential links between large-scale drivers and the variability of precipitation extremes at the regional and local scales. Caesar et al. (2011) pointed out the possibility of a positive correlation between several extreme precipitation indices with a La Niña-like sea surface temperature pattern over the Indo-Pacific region. Ho et al. (2011) also pertained to potential changes in heavy rainfall over Vietnam in relationship with the sea surface temperature patterns over the South China Sea. This study further indicated that the annual peaks of several extreme rainfall indices, characterizing intensity, frequency, and wet duration, usually coincided with negative IOD and La Niña events, while dry extremes were usually found during El Niño events. Further studies should pay more attention to large-scale controls on non-monotonic trend components in extreme precipitation indices.

Conclusions
It is always a good practice to perform a comprehensive evaluation of systematic changes in weather and climate extremes, which is conducive to local decision-makers to put forward more informed decisions on climate change adaptation and mitigation strategies. Thus, the present work aimed to assess and estimate observed changes in heavy rainfall events over the Central Highlands of Vietnam during recent decades. On this basis, ten ETCCDI extreme rainfall indices (i.e., Rx1day, Rx5day, SDII, R95p, R99p, R10mm, R20mm, R50mm, CDD, and CWD) were subjected to temporal trend analysis by employing an innovative trend analysis methodology. The ITA method can identify non-monotonic trend components in a given hydro-meteorological time series without considering any restrictive assumptions. Such an easy-to-apply but productive technique can be utilized to explore trend possibilities in a wide range of science and engineering fields.
As expected, observed trends in extreme rainfall indices are less spatially consistent across the Central Highlands irrespective of overall or partial trends. Moreover, most analyzed extreme indices exhibited non-monotonic forms, in which low-, medium-, and high-value subgroups exposed opposite trend directions or the same tendency with very different magnitude. In particular, the high-value clusters were characterized mainly by significant increases in most analyzed extreme indices such as Rx5day, SDII, R95p, R99p, R50mm, and CDD. It is worth pointing out that several stations also exposed significant increasing trends within all low-, medium-, and high-value subgroups of Rx5day, R95p, R99p, R20mm, and R50mm (i.e., Kon Tum), SDIII (i.e., Bao Loc), and CDD (i.e., Ayun Pa). On the whole, these findings imply significant increases in the frequency and intensity of heavy rainfall events over the Central Highlands of Vietnam during recent decades, which should be taken into account when devising any local adaptation actions to climate-related risks in the study area. Asterisks "*" indicate that the correlation is statistically significant at the 95% confidence level Spearman's rank correlation coefficient (rho)