Spatiotemporal Variability of Precipitation and Teleconnections in Mekong Delta (Vietnam)

1. Introduction

Precipitation variability is a critical determinant of water resource availability, agricultural productivity, and socioeconomic stability in tropical monsoon regions. In Southeast Asia, interannual and interdecadal rainfall fluctuations are primarily influenced by large-scale ocean–atmospheric teleconnection patterns. Among these, the El Niño–Southern Oscillation (ENSO), Indian Ocean Dipole (IOD), and Pacific Decadal Oscillation (PDO) are the principal drivers [1,2,3]. A comprehensive understanding of how these climate modes, both individually and collectively, modulate regional precipitation is essential for effective climate adaptation and water management planning in the region.

The VMD, situated at the downstream end of the Mekong River system and spanning approximately 39,000 km2 in southern Vietnam, is one of the most hydrologically complex and economically significant regions in Southeast Asia. As Vietnam's principal agricultural area, the VMD accounts for more than 50% of the national rice production and nearly 70% of the fruit and aquaculture yields. Despite its productivity, the region faces significant exposure to climate extremes, including seasonal flooding from upstream Mekong inflows and severe dry-season droughts with saline intrusion along the coastline. These contrasting hydrological regimes are strongly influenced by interannual climate variability; however, the underlying teleconnection mechanisms remain insufficiently understood at the subregional scale.

ENSO, defined by anomalous sea surface temperature (SST) warming (El Niño) or cooling (La Niña) in the central and eastern equatorial Pacific, significantly influences global precipitation patterns by altering the Walker Circulation and Hadley Cell [4,5]. In Southeast Asia, El Niño events are typically linked to reduced monsoon rainfall and increased drought risk, whereas La Niña events often result in enhanced precipitation and potential flooding [6]. In Vietnam, ENSO-driven precipitation anomalies have been observed across multiple river basins, with El Niño years showing rainfall deficits of 10–30% in the central and southern regions [7].

The IOD, quantified by the DMI as the SST gradient between the western (50–70°E) and eastern (90–110°E) tropical Indian Ocean [3], represents an independent but often ENSO-coupled source of climate variability. Positive IOD events, characterized by anomalously warm western and cool eastern Indian Ocean SSTs, tend to weaken cross-equatorial moisture transport toward Southeast Asia, thereby reducing precipitation over the Indo-Chinese Peninsula [8,9]. Critically, co-occurring positive ENSO and positive IOD events, such as those observed in 1997–1998, 2015, and 2019, have been shown to produce compound drying effects that are substantially more severe than either mode in isolation [10], underscoring the importance of multi-index analyses.

The PDO, defined as the leading mode of North Pacific sea surface temperature (SST) variability on decadal timescales [11], does not directly influence tropical precipitation patterns. Instead, it serves as a low-frequency modulator of ENSO amplitude and frequency [12]. Positive PDO phases are associated with more frequent and intense El Niño events, which amplify the drying effects on Southeast Asian monsoon rainfall. In contrast, negative PDO phases may mitigate the El Niño-induced rainfall deficits. This decadal background state is therefore essential for interpreting the non-stationarity of ENSO teleconnections observed across different decades [13,14].

Although research on climate teleconnections in Southeast Asia is expanding, significant knowledge gaps remain regarding VMD. First, most existing studies have examined ENSO effects on national or regional precipitation using spatially aggregated datasets. This approach obscures the distinct hydrological gradient between the upper delta, which is primarily influenced by Mekong River inflows and continental rainfall, and the coastal margins, which are affected by tidal forcing, saline intrusion, and oceanic moisture from the East Sea and Gulf of Thailand [15,16]. Second, few studies have jointly assessed the contributions of all three major climate modes (ENSO, IOD, and PDO) to precipitation variability. Most investigations focus exclusively on ENSO [17,18], neglecting the compounding influence of the IOD and the modulating effect of the PDO. This single-index approach is inadequate, given the documented interactions among these modes and their nonlinear impacts on regional hydrology. Third, the time-varying nature of teleconnection strengths, such as the possibility that ENSO's influence on VMD rainfall may intensifies during certain decades or at specific frequency bands, has not been systematically investigated using time-frequency methods such as WTC. While standard Pearson correlation analysis is useful for quantifying mean linear relationships and optimal lag structures, it cannot capture the non-stationary and multi-scale characteristics of ocean–atmosphere–precipitation linkages [19,20].

This study sought to deliver a comprehensive, multi-index, and multi-scale analysis of the teleconnection between large-scale climate variability and precipitation anomalies (PPTA) in the Vietnamese Mekong Delta from 1981 to 2025. The specific objectives were as follows:

(1) To characterize the spatiotemporal variability of monthly precipitation across ten meteorological stations representing three hydrological zones of the Vietnamese Mekong Delta: the upper delta, eastern coastal margin, and southern and western coastal margins.

(2) To quantify the lag-correlation structure between PPTA and the three climate indices (Niño 3.4, DMI, and PDO) at both the station and cluster levels, and to identify optimal lag periods and spatial gradients in teleconnection strength.

(3) To analyze the time-varying, frequency-specific coherence between climate indices and PPTA using Wavelet Transform Coherence (WTC), in order to reveal non-stationary teleconnection patterns across interannual (2–5 year) and interdecadal (8–16 year) bands.

(4) To elucidate the differential sensitivity of hydrological clusters within the Vietnamese Mekong Delta to each climate mode and discuss the implications for drought early warning, water resource management, and climate-smart agriculture in the region.

This study makes several novel contributions. This is the first study to examine the teleconnection signatures of ENSO, IOD, and PDO on sub-regional precipitation variability across the VMD using a 45-year high-resolution station dataset. By integrating Pearson lag-correlation analysis with wavelet transform coherence across three distinct hydrological clusters, this study provides a spatially differentiated understanding of climate forcing that informs decision-making for regional water managers and agricultural planners. Furthermore, conceptualizing the PDO as a low-frequency background modulator, rather than a direct forcing agent, advances the framework for interpreting decadal shifts in ENSO teleconnection strength over the delta. These findings are particularly relevant as climate change is projected to intensify both El Niño and La Niña events, increasing the risks to the VMD's water security and food production systems.

2. Materials and Methods

2.1. Study Area

The VMD is located in the southernmost part of Vietnam, bounded by latitudes 8.5°N to 11.0°N and longitudes 104.5°E to 106.8°E (Figure 1). This delta is formed by the terminal distributary network of the Mekong River, which discharges into the East Sea through the Tien and Hau Rivers, and into the Gulf of Thailand via the Ca Mau Peninsula. The region's flat and low-lying topography, with elevations rarely exceeding 5 m above sea level, makes it highly susceptible to both riverine flooding and tidal inundation from the south and west [16,21].

The VMD experiences a tropical monsoon climate characterized by two main seasons: a wet season from May to November, driven by the southwest monsoon and accounting for approximately 85–90 percent of annual rainfall, and a pronounced dry season from December to April, marked by minimal precipitation and increased saline intrusion along the coastal fringe [15,22]. The mean annual precipitation across the delta varies significantly, ranging from approximately 1,370 mm at inland stations such as Cao Lanh to over 1,820 mm at the southernmost tip, Ca Mau. This variation reflects the influence of orographic effects, proximity to coastal moisture sources, and exposure to cyclonic systems in the Gulf of Thailand.

Owing to the pronounced spatial heterogeneity in hydrological regimes and climate sensitivity, the ten selected meteorological stations were classified into three ecologically distinct hydrological clusters:

• Upper Delta Cluster (Chau Doc, Cao Lanh, Can Tho): Located in the northwestern and central delta, these stations receive Mekong River floodwaters directly from Cambodia. They are characterized by alternating flood and drought regimes, which are primarily driven by upstream discharge. The mean annual precipitation ranges from 1,373 to 1,535 mm.
• Eastern coastal cluster (My Tho, Cang Long, Ba Tri): Situated along the northeastern delta margin facing the East Sea, these stations are subject to strong tidal influence through the Tien and Ham Luong River mouths. They exhibit increased rainfall variability associated with northeast monsoon interaction. The mean annual precipitation ranges from 1,406 to 1,587 mm.
• Southern–Western coastal cluster (Soc Trang, Bac Lieu, Ca Mau, and Rach Gia): encompassing the southernmost and western coastal margins bordering the Gulf of Thailand, these stations recorded the highest annual rainfall totals (1,636–1,823 mm). They also demonstrated the greatest vulnerability to dry-season droughts, large-scale shrimp aquaculture exposure, and saline intrusion. Ca Mau, the southernmost station, recorded the highest drought frequency in the delta [23].

2.2. Precipitation Data

Monthly precipitation totals (mm/month) for the ten stations were sourced from the Southern Regional Hydrometeorological Center (SRHMC) of the Vietnam Meteorological and Hydrological Administration (VNMHA) for the period January 1981 to December 2025, resulting in a continuous record of 540 monthly observations per station. Daily station records underwent quality control procedures, including range tests that excluded values exceeding five standard deviations from the long-term monthly mean and internal consistency checks prior to aggregation into monthly totals. A monthly record was designated as missing if more than ten days within that month lacked valid observations. The resulting monthly precipitation dataset demonstrated mean annual totals consistent with previously published climatological analyses for the [15,23].

2.3. Climate Indices

Three large-scale climate indices were employed to represent the principal modes of interannual and interdecadal ocean–atmosphere variability relevant to precipitation in Southeast Asia (Table 1).

Niño 3.4 index: The Niño 3.4 SST anomaly index, defined as the area-averaged SST departure from the 1991–2020 climatological baseline over the region 5°N–5°S, 170°W–120°W, was obtained from the NOAA Extended Reconstructed Sea Surface Temperature dataset version 5 (ERSST v5). This index is the most widely adopted indicator of ENSO phase and intensity, with values exceeding +0.5°C and −0.5°C sustained over five consecutive overlapping three-month periods defining El Niño and La Niña conditions, respectively [24,25].

Dipole Mode Index (DMI): The DMI quantifies the IOD as the difference in area-averaged SST anomalies between the western Indian Ocean (50–70°E, 10°S–10°N) and the eastern Indian Ocean (90–110°E, 10°S–0°N) [3]. The monthly DMI values were sourced from the Hadley Centre Sea Ice and Sea Surface Temperature (HadISST) dataset, accessed via the NOAA Physical Sciences Laboratory. Positive DMI events (DMI > +0.4°C) are associated with anomalously warm western and cool eastern Indian Ocean conditions, which suppress moisture flux toward Southeast Asia and reduce precipitation over the VMD region.

Pacific Decadal Oscillation (PDO): The PDO index, representing the leading empirical orthogonal function of monthly North Pacific SST anomalies north of 20°N [11], was obtained from the NOAA PSL using the ERSST v5 dataset. The PDO captures decadal-scale SST variability (typical period 15–30 years) and serves as a modulator of ENSO teleconnection strength. In its positive phase, the PDO amplifies El Niño-like warming in the tropical Pacific, potentially intensifying ENSO-driven precipitation deficits in Southeast Asia [12,14].

All three climate index time series covered the same 1981–2025 period as the precipitation records, providing 540 monthly observations. Prior to analysis, the following preprocessing steps were applied consistently across all datasets: (i) standardization of temporal alignment to ensure that all series shared the same monthly time step referenced to the first day of each month; (ii) replacement of documented fill values (e.g., −9999) with not-a-number (NA) markers; (iii) linear interpolation of isolated missing values, with monthly climatological means substituted where gaps exceeded two consecutive months; and (iv) verification that the range of each index was physically plausible (|Niño 3.4| < 4°C; |DMI| < 2°C; |PDO| < 5). The integrated master dataset, comprising 540 rows and 13 variables (Year, Month, ten station precipitation columns, and three climate index columns), was used as the sole input for all statistical and spectral analyses.

2.4. Methods

The analytical framework of this study comprised three sequential components: (i) descriptive characterization of precipitation variability and PPTA computation; (ii) Pearson lag-correlation analysis to quantify the linear, time-averaged relationship between climate indices and PPTA at the station and cluster levels; and (iii) Wavelet Transform Coherence (WTC) analysis to reveal the time-varying and frequency-specific nature of teleconnection signals. A summary of the methods, their purposes, and key references are provided in Table 2.

2.4.1. Precipitation Anomaly Computation

Monthly precipitation anomalies were expressed as the Precipitation Anomaly Percentage (PPTA, %) relative to the long-term monthly climatological mean computed over the full study period (1981–2025). The precipitation anomaly percentage (PPTA, %) was computed for each station and month as follows:

$${PPTA}_{(i,m)}={\displaystyle \frac{\left[{PPT}_{(i,m)}-{PP\overline{T}}_{\left(m\right)}\right]}{{PP\overline{T}}_{\left(m\right)}}}\times 100$$

(1)

where PPT_(i,m) is the observed precipitation at station i in calendar month m, and PPT̄_(m) is the long-term monthly mean computed over the full baseline period of 1981–2025. A regional PPTA index was derived as the arithmetic mean of PPTA across all ten stations for each month, providing a spatially integrated measure of delta-wide precipitation anomalies.

$${PPTA\_regional}_{(m,y)}={\displaystyle \frac{1}{N}}\sum _i{PPTA}_{(i, m,y)}$$

(2)

where N = 10 is the number of stations, and Annual precipitation anomalies were derived by summing the monthly PPTA values within each calendar year. Positive PPTA indicates above-normal rainfall relative to the climatological baseline, whereas negative PPTA indicates below-normal conditions associated with potential agricultural drought.

2.4.2. Pearson Lag-Correlation Analysis

1.

Correlation coefficient

The Pearson product-moment correlation coefficient (r) was used to quantify the linear relationship between each climate index and PPTA series. Given two stationary time series X = {x₁, x₂, …, xₙ} (climate index) and Y = {y₁, y₂, …, yₙ} (PPTA), the Pearson correlation is defined as

$$r(X,Y)={\displaystyle \frac{\sum \left[\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)\right]}{\sqrt{{\left(x_i-\overline{x}\right)}^2{\left(y_i-\overline{y}\right)}^2}}}$$

(3)

where x̄ and ȳ denote the sample means of X and Y respectively. The value of r ranges from −1 (perfect negative linear relationship) to +1 (perfect positive linear relationship), with values near zero indicating no linear association [27,28]. Statistical significance was evaluated using a two-tailed t-test with (n − 2) degrees of freedom at the 5% (p < 0.05) and 1% (p < 0.01) significance levels.

2)

Lag-correlation structure

To account for the delayed response of regional precipitation to large-scale SST anomalies, reflecting the time required for oceanic and atmospheric teleconnection signals to propagate from the Pacific and Indian Oceans to mainland Southeast Asia, lag-correlation analysis was performed for lags l = 0, 1, 2, …, 12 months. At each lag, the climate index series X was shifted forward by l months relative to the PPTA series Y, and the Pearson correlation was computed between the truncated series as follows:

$$r\left(l\right)=Pearson\left[X\left(t\right),Y(t+1)\right] for l=0, 1, \dots , 12$$

(4)

The optimal lag l* is defined as the lag at which the absolute correlation |r(l)| is maximized and statistically significant (p < 0.05). This procedure was applied independently for each combination of climate index (Niño 3.4, DMI, PDO) and PPTA series (10 individual stations + 1 regional mean), yielding a comprehensive lag correlation matrix. The results were additionally aggregated by hydrological cluster to identify spatial gradients in teleconnection strength and lag structure across the three ecological zones of the VMD.

3)

Severity-stratified correlation

Following the approach of [26], Pearson correlations were computed separately for distinct ENSO and IOD severity classes to assess whether teleconnection strength varied with the intensity of the forcing event. ENSO phases were classified using the Niño 3.4 index relative to the ±0.5°C threshold: weak (±0.5 to ±1.0°C), moderate (±1.0 to ±1.5°C), and strong (|Niño 3.4| > 1.5°C). IOD phases were similarly classified using the DMI: weak (±0.2 to ±0.4°C), moderate (±0.4 to ±0.8°C), and strong (|DMI| > 0.8°C) [3]. For each severity class, the correlation between the respective index and regional PPTA was computed using all months that satisfied the classification criterion.

2.4.3. Wavelet Transform Coherence Analysis

4)

Continuous wavelet transform

Wavelet analysis simultaneously decomposes a time series in the time and frequency domains, enabling the identification of transient periodicities and their temporal evolution [20,29]. The Continuous Wavelet Transform (CWT) of a discrete time series xₙ (n = 1, 2, …, N) with a uniform time step δt is defined as the convolution of xₙ with a scaled and translated version of a mother wavelet ψ₀(η):

$$W_n^x\left(S\right)= \sqrt{\delta \left({\displaystyle \frac{t}{s}}\right){\sum }_n^{\text{'}}{x}_n^{,}{\psi }^{*}\left[(n^{,}-n)\delta \left({\displaystyle \frac{t}{s}}\right)\right]}$$

(5)

where s is the wavelet scale, n is the localized time index, and ψ* denotes the complex conjugate of the mother wavelet, respectively. The Morlet wavelet was selected as the mother wavelet for this study because it provides an optimal balance between time localization and frequency resolution. The Morlet wavelet is defined as

$${\psi }^0\left(\eta \right)={\pi }^{\left(-{\displaystyle \frac{1}{4}}\right)}\exp \left(i{\varpi }_0\eta \right)exp\left({-\eta }^2/2\right)$$

(6)

where ω₀ = 6 is the non-dimensional angular frequency that satisfies the admissibility condition [29]. At ω₀ = 6, the Morlet wavelet has a central frequency of approximately one cycle per unit scale, and the wavelet scale s is approximately equal to the Fourier period, facilitating the physical interpretation of the results.

2.

Wavelet transform coherence

To examine the co-variability between two time series, a climate index X and a PPTA series Y, in both time and frequency domains simultaneously, Wavelet Transform Coherence (WTC) was employed [19,30]. The WTC is defined as the ratio of the smoothed cross-wavelet spectrum to the product of the smoothed individual power spectra:

$$R_n^2\left(s\right)={\displaystyle \frac{{\left|S\left(s^{-1}{W}_n^x\right(s\left){W_n^y}^{*}\right(s\left)\right)\right|}^2}{\left[S\left(s^{-1}{\left|W_n^x\left(s\right)\right|}^2\right).S\left(s^{-1}{\left|W_n^y\left(s\right)\right|}^2\right)\right]}}$$

(7)

where S(·) is a smoothing operator applied in both the time and scale dimensions, Wˣₙ(s) and Wʸₙ(s) are the CWT coefficients of X and Y, and the asterisk (*) denotes the complex conjugate. The WTC value R²ₙ(s) ranges between 0 and 1, analogous to a squared correlation coefficient in the time-frequency domain: values close to 1 indicate strong co-variability between the two series at scale s and time n, while values near 0 indicate weak or absent coherence [31].

The smoothing operator S(·) in Eq. (7) decomposes into time-domain smoothing Sₜᵢₘₑ and scale-domain smoothing Sₛₒₐₗₑ, implemented as

$$S_{time}{\left.\left(W\right)\right|}_s=\left[W_n(s)\times c_1\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \frac{{-t}^2}{{2s}^2}}\right)\right] (smoothing in time)$$

(8a)

$$S_{scale}{\left.\left(W\right)\right|}_n=\left[W_n\left(s\right)\times c_2\mathsf{\Pi }\left(0.6s\right)\right] (smoothing in scale)$$

(8a)

where c₁ and c₂ are normalization constants, and Π denotes a rectangular window of width 0.6s in scale space, following the implementation described by [30] and [19].

3.

Phase analysis

The relative phase angle between the two series was extracted from the complex cross-wavelet spectrum and overlaid on the WTC as a directional arrow. Phase relationships are interpreted as follows: rightward-pointing arrows (→) indicate in-phase co-variability (positive correlation), where both series increase or decrease simultaneously; leftward-pointing arrows (←) indicate anti-phase (negative correlation), where one series increases as the other decreases; upward arrows (↑) indicate that X leads Y by one-quarter cycle; and downward arrows (↓) indicate that Y leads X by one-quarter cycle [19]. Phase arrows were plotted only where R²ₙ(s) > 0.7 to highlight the regions of strong coherence.

4.

Statistical significance testing

Because the theoretical distribution of the WTC under the null hypothesis of no coherence cannot be derived analytically, statistical significance was assessed using a Monte Carlo simulation approach [19,20]. For each pair of time series, 300 surrogate pairs were generated by phase randomization of the original series, preserving their power spectra but destroying any temporal co-variability structure. The WTC was computed for each surrogate pair, and the 95th percentile of the resulting WTC distribution was adopted as the significance threshold at each scale and time position. Regions where the observed R²ₙ(s) exceeds this threshold are enclosed by white significance contours in the WTC plots and are interpreted as representing statistically significant coherence at the 95% confidence level.

Analyses were performed for three representative stations: Chau Doc (upper delta), Ba Tri (eastern coastal), and Ca Mau (southern–western coastal) as proxies for each hydrological cluster, yielding nine WTC spectra (3 stations × 3 indices). The cone of influence (COI), within which edge effects may distort the wavelet power, was computed following [20] and delineated on all WTC plots. Periodicities outside the COI were excluded from the interpretation. All wavelet computations were performed using the biwavelet package (v0.20.21) in R [32].

2.5. Analytical Framework and Software

The complete analytical pipeline, from raw data ingestion through PPTA computation, Pearson lag-correlation, and WTC analysis, was implemented in R (version 4.3.x) using the following packages: dplyr and tidyr for data manipulation [33]; ggplot2 for visualization, zoo for time series interpolation; and biwavelet for CWT and WTC computation.

A summary of the complete methodological workflow is as follows: (1) quality-controlled monthly precipitation data from 10 VMD stations and three climate indices (Niño 3.4, DMI, PDO) spanning 1981–2025 were assembled into an integrated master dataset; (2) PPTA was computed for each station and aggregated into cluster and regional indices; (3) Pearson lag-correlation was performed for lags 0–12 months between each index and PPTA pair, with results stratified by ENSO/IOD severity class; (4) WTC was applied to representative station–index pairs, with Monte Carlo significance testing and phase arrow visualization; and (5) results were synthesized across methods to identify convergent teleconnection signals and their spatial differentiation across the VMD's hydrological gradient. A schematic overview of the proposed methodology is presented in Figure 2.

3. Results

3.1. Spatiotemporal Variability of Precipitation

3.1.1. Climatological Distribution

The mean annual precipitation at the ten VMD stations during 1981–2025 ranged from 1,373 mm at Cao Lanh in the upper delta to 1,823 mm at Ca Mau on the southern coast, with a regional mean of approximately 1,578 mm per year (Table 3). A pronounced north–south and inland–coastal gradient was observed. The southern–western coastal cluster (Soc Trang, Bac Lieu, Ca Mau, and Rach Gia) consistently recorded higher annual precipitation totals (1,636–1,823 mm) than the upper delta cluster (1,373–1,535 mm), which reflects the combined effects of Gulf of Thailand moisture advection and orographic enhancement along the Ca Mau Peninsula.

The monthly precipitation distribution demonstrated a distinctly unimodal wet season. The mean monthly rainfall remains low from December to February (8.5–19.9 mm/month), increases markedly beginning in April with the onset of the southwest monsoon, and reaches its maximum between June and October (223–240 mm/month) before declining rapidly in November. The wet season (May–November) contributes 87–91% of the annual rainfall at all stations. The highest wet-season proportion was recorded at Ca Mau (91.0%), whereas the lowest was recorded at Chau Doc (88.2%), reflecting the greater dry-season precipitation at the upper delta stations influenced by residual Mekong River moisture.

Intra-annual variability, measured by the coefficient of variation (CV) of monthly rainfall, was inversely correlated with mean precipitation. Southern coastal stations, such as Ca Mau and Bac Lieu, displayed lower CV values (68.9–72.1%) than upper delta stations, including Chau Doc (78.5%) and Cao Lanh (82.3%). This pattern indicates that stations with higher rainfall generally experience more consistent precipitation distribution throughout the year.

3.1.2. Interannual Variability and PPTA

The annual PPTA series for the regional mean (Figure 3) revealed pronounced interannual variability over the 1981–2025 period.

The long-term mean annual precipitation is 1,578 mm per year. Over the study period, 44% of the years experienced above-normal annual rainfall (PPTA > 0%), whereas 56% exhibited below-normal conditions. The most extreme above-normal year occurred in 2020 (PPTA = +70.5%), whereas the most extreme below-normal year was 1982 (PPTA = −73.0%), coinciding with the strong 1982–1983 El Niño event (annual mean Niño 3.4 = +0.83°C).

Persistent multi-year periods of below-normal precipitation occurred during 1982–1991 and 2001–2015, whereas extended wet periods were observed during 1992–2000 and 2016–2025. The transition from below-normal to above-normal regimes around 2016 is particularly notable, coinciding with a shift in the PDO phase from positive to negative. This pattern aligns with the established decadal modulation of Southeast Asian monsoon rainfall [12,14]. Years associated with major El Niño events (1982, 1987, 1997, and 2015) generally exhibited negative PPTA, whereas La Niña years (1996, 1999, 2021, and 2022) were predominantly linked to positive anomalies.

3.2. Pearson Lag-Correlation Analysis

3.2.1. Precipitation Anomaly Computation

Figure 4 presents the lagged Pearson correlation coefficients between large-scale climate indices (Niño 3.4, DMI, and PDO) and regional precipitation anomalies (PPTA) in the Mekong Delta over a 0–12 month lag period. The red dashed lines denote the 95% confidence interval (p < 0.05), which highlights statistically significant relationships. These findings demonstrate a pronounced lag-dependent teleconnection structure with distinct temporal responses associated with each climate mode.

Niño 3.4: a strong and consistently negative correlation exists between Niño 3.4 and regional PPTA across lags of 0 to 7 months, with all correlations significant at the 1% level (p < 0.01). The correlation reached its maximum magnitude at a lag of two months (r = −0.304, p < 0.001), suggesting that ENSO sea surface temperature anomalies in the central Pacific precede VMD precipitation anomalies by approximately two months. This lag corresponds to the propagation time required for atmospheric Rossby wave adjustments from the Pacific to the Indo-Chinese Peninsula, mediated by changes in the Walker Circulation and the South Asian monsoon trough [25]. The negative correlation indicates that El Niño events (positive Niño 3.4) are associated with reduced VMD rainfall, whereas La Niña events (negative Niño 3.4) are associated with increased rainfall.

DMI: The relationship between DMI and PPTA is notably weaker and exhibits a distinct temporal pattern. Correlations at lags of 0 to 7 months were negligible and lacked statistical significance. A significant positive correlation was observed at lags of 8 to 12 months, with the highest values at lag 11 (r = 0.186, p < 0.001) and lag 12 (r = 0.182, p < 0.001). This positive correlation at longer lags may indicate a delayed response of Indian Ocean sea surface temperature gradients to ENSO-driven thermocline adjustments, which, in turn, influence the timing of the following year's monsoon onset over Southeast Asia [8,10]. The relatively weak DMI signal at the regional level likely results from partial cancellation among individual stations within the three distinct clusters, each of which is variably influenced by Indian Ocean moisture pathways.

PDO: the PDO demonstrated a persistent and statistically significant negative correlation with regional precipitation anomalies (PPTA) across all 13 lag values (lag 0–12), with correlation coefficients ranging from −0.094 to −0.253 (all p < 0.05). This broad, lag-insensitive negative association aligns with the PDO's function as a low-frequency modulator of the winter climate. Positive PDO phases are associated with increased El Niño frequency and intensity, which suppresses VMD precipitation by modulating the background Pacific sea surface temperature (SST) state rather than through a sharply defined teleconnection pathway. The negative correlation was most pronounced at lags of 2–5 months (r = −0.25), indicating a potential contribution from seasonal-scale Pacific–Indian Ocean interactions.

3.2.2. Station-Level Analysis and Spatial Gradient

Table 4 presents the optimal lag and corresponding Pearson correlations for each station–index combination. The spatial pattern was broadly consistent across stations but revealed systematic differences across hydrological clusters.

Niño 3.4: The optimal lag for Niño 3.4 was two months at the upper delta and certain eastern coastal stations (Chau Doc, Cao Lanh, Can Tho, My Tho, and Rach Gia), and three months at the eastern and southern coastal stations (Cang Long, Ba Tri, Soc Trang, Bac Lieu, and Ca Mau). This 1-month difference in the optimal lag indicates that the ENSO teleconnection signal reaches the upper delta more rapidly, possibly through direct continental moisture advection, than the southern coastal margins, which experience additional oceanic modulation from the Gulf of Thailand. Correlation magnitudes were highest at Ca Mau (r = −0.302) and lowest at Chau Doc (r = −0.271), reflecting a somewhat stronger ENSO sensitivity at southern coastal stations, consistent with their increased exposure to maritime moisture sources.

DMI: The optimal DMI lag was 11 months at the upper delta and most eastern coastal stations, and 12 months at the three southernmost stations (Soc Trang, Bac Lieu, and Ca Mau), indicating a slightly longer Indian Ocean teleconnection pathway to the southern coast. Although this pattern is consistent, the DMI correlation magnitudes are substantially weaker than those of Niño 3.4 across all stations (r = 0.140–0.195). While statistically significant at the 1% level, these correlations account for a considerably smaller proportion of the PPTA variance. These findings indicate that the IOD exerts a secondary, yet non-negligible, modulating influence on VMD precipitation, primarily through lagged SST-monsoon interactions in the Arabian Sea.

PDO: The optimal lag for the PDO ranged from 2 months in the upper delta and most stations to 4 months at Ba Tri and Ca Mau, with correlation magnitudes between −0.207 and −0.263. The longer optimal lag observed at the southern coastal stations suggests varying exposure to the Pacific and Gulf of Thailand moisture pathways. The PDO signal was strongest in the upper delta (Cao Lanh: r = −0.244) and eastern coastal (Cang Long: r = −0.263) clusters, and weakest at the westernmost station, Rach Gia (r = −0.207).

Figure 5 displays scatter plots of regional precipitation anomalies (PPTA) in relation to the three climate indices at their respective optimal lag times, as determined by lag-correlation analysis: Niño 3.4 (lag of 2 months), DMI (lag of 11 months), and PDO (lag of 2 months). Linear regression lines and Pearson correlation coefficients were included to quantify the strength and direction of these associations. These scatter plots corroborate the lag-correlation analysis by depicting the direction, magnitude, and dispersion of teleconnection signals at their most influential timescales. ENSO was identified as the primary driver at short lags, PDO as a secondary but persistent modulator, and DMI as a delayed and weaker influence. The considerable scatter in all panels indicates that precipitation variability in the Mekong Delta results from both large-scale climate forcing and local or regional processes.

3.3. Wavelet Transform Coherence Analysis

3.3.1. Niño 3.4 – PPTA Coherence

Figure 6, Figure 7 and Figure 8 illustrate the WTC between Niño 3.4 and PPTA at the three representative stations. At all stations, significant coherence (R2 > 0.5) was concentrated within the 2–5 year period band, aligning with the established 2–7 year recurrence interval of ENSO events [2,20]. Phase arrows in significant regions were predominantly leftward-pointing, corroborating the anti-phase relationship identified by Pearson’s analysis: El Niño conditions suppressed VMD precipitation, whereas La Niña conditions enhanced it.

Coherence is temporally intermittent rather than continuous, with the most prominent significant regions observed during 1985–1995 (encompassing the 1987–1988 El Niño and the strong 1988–1989 La Niña), 1997–2000 (the major 1997–1998 El Niño–La Niña transition), and 2014–2018 (the 2015–2016 strong El Niño) periods. This episodic pattern highlights the non-stationary nature of the ENSO–VMD precipitation teleconnection, which is not fully captured by time-averaged correlations.

Spatial analysis indicates that Niño 3.4 WTC coherence is strongest and most persistent at Ca Mau (southern coast), where significant coherence at 2–4 year periods extended from approximately 1985 to 2020. At Chau Doc (upper delta), coherence was weaker and more episodic, reflecting the buffering influence of upstream Mekong River inflows and continental precipitation processes on the local ENSO signal. Ba Tri (eastern coastal) demonstrated intermediate coherence strength, with notable power at 2–4 year periods during 1990–2005.

3.3.2. DMI – PPTA Coherence

The WTC analysis between the DMI and PPTA (Figure 9, Figure 10 and Figure 11) demonstrates a coherence structure that differs markedly from that of Niño 3.4. The significant coherence regions were generally weaker and more fragmented, which aligns with the observed weaker Pearson correlations. Nevertheless, distinct coherence patches were observed in the 1–2 year period band from 2000 to 2010 and in the 4–6 year band from 2008 to 2018. These patches are characterized predominantly by rightward-pointing arrows, indicating an in-phase or leading-phase relationship between the DMI and PPTA at these temporal scales.

The positive in-phase relationship at extended lags (11–12 months, as identified in the Pearson analysis) is evident in the WTC as a slight phase lead structure within the annual period band. Upward and rightward arrows indicate that DMI variations from the previous year were positively associated with PPTA in the subsequent year. This observation aligns with the proposed mechanism in which negative IOD events, characterized by cooler eastern Indian Ocean sea surface temperatures, enhance moisture flux to Southeast Asia during the following wet season [3,8]. The IOD signal was most pronounced at southern and western coastal stations, such as Ca Mau and Rach Gia, which are more directly influenced by moisture pathways from the Gulf of Thailand that are affected by Indian Ocean sea surface temperature gradients.

3.3.3. PDO – PPTA Coherence

The PDO–PPTA wavelet transform coherence (WTC) (Figure 12, Figure 13 and Figure 14) demonstrates significant coherence primarily within the 8–16 year period band, aligning with the Pacific Decadal Oscillation's (PDO) characteristic decadal-to-interdecadal variability timescale [11,12]. Notably, significant anti-phase coherence (indicated by leftward arrows) at periods of 10–14 years is evident during approximately 1985–2010 at all three representative stations. This pattern reflects the transition from the negative PDO phase (1947–1976) to the positive PDO phase (1977–2000), followed by a partial return to the negative phase. Such decadal-scale coherence accounts for the broad, lag-insensitive negative Pearson correlation identified in Section 3.2.1 and substantiates the interpretation of the PDO as a background modulator that amplifies or attenuates the interannual ENSO signal in the VMD.

At shorter periods (2–4 years), the PDO–PPTA WTC exhibits intermittent yet statistically significant coherence during 1997–2005 and 2013–2020 periods. This likely reflects the overlap between PDO-modulated Pacific sea surface temperature (SST) and concurrent ENSO events. The mixed-phase structure at these shorter periods suggests that the PDO's influence on sub-decadal precipitation variability in the VMD is mediated through ENSO rather than via a direct atmospheric pathway [14].

3.3.3. Summary of Teleconnection Signals

The integration of Pearson lag-correlation and WTC analyses revealed three principal teleconnection signals influencing VMD precipitation variability. First, a dominant interannual ENSO signal operated primarily at 2–5 year periodicities with a 2–3 month lag, exhibiting the greatest strength at the southern coastal cluster. Second, a secondary IOD signal displays a long lagged response (11–12 months) and weaker coherence, with effects most pronounced at the southern and western coastal margins. Third, a background PDO modulation occurs at decadal (8–16 year) scales, remains broadly consistent across all three clusters, and manifests as persistent low-frequency damping or amplification of ENSO teleconnection. The combined multi-index analysis demonstrated that no single climate mode fully explains the VMD precipitation variability. Compound events, particularly those with simultaneous El Niño and positive IOD phases (1997, 2015, 2019), resulted in the most severe and spatially coherent precipitation deficits across the entire delta.

4. Discussion

The results presented in Section 4 provide a comprehensive multi-scale characterization of the teleconnection between large-scale climate variability and precipitation anomalies across the VMD over the 1981–2025 period. This section contextualizes the findings within the broader literature on Southeast Asian climate–precipitation linkages, discusses the mechanistic pathways underlying the observed teleconnections, examines the implications of compound climate forcing and PDO modulation, interprets the spatial gradients across hydrological clusters, and identifies limitations and directions for future research.

4.1. ENSO Teleconnection: Mechanisms, Lag Structure, and Comparison with Regional Studies

The principal finding of this research is a statistically significant negative correlation between the Niño 3.4 index and VMD PPTA, with the strongest relationship observed at a 2-month lag (r = −0.304, p < 0.001) at the regional scale and −0.302 at Ca Mau with a 3-month lag. This outcome aligns with previously documented ENSO–precipitation teleconnections across mainland Southeast Asia [7,13,18], and substantiates that El Niño suppresses, while La Niña enhances, VMD precipitation by altering the Walker Circulation and the intensity of the South Asian monsoon trough [25].

The 2–3 month optimal lag identified in this study is shorter than the 4–5 month lag reported by [26] for northeast Thailand and the 4-month lag documented by [6]. This difference is physically significant because the VMD is situated approximately 800–1,200 km south of northeastern Thailand, placing it closer to the equatorial Pacific and the primary Walker Circulation anomalies. A shorter lag indicates that ENSO-driven atmospheric teleconnections, specifically anomalous subsidence over Southeast Asia resulting from El Niño's eastward-shifted convective heating, reach the VMD region more rapidly via equatorial Kelvin and Rossby wave propagation [34]. Furthermore, the VMD's direct exposure to the South China Sea and Gulf of Thailand, which act as conduits for ENSO-modulated moisture flux from the western Pacific, may further accelerate atmospheric response.

The Niño 3.4 index accounted for approximately 9.2% of the total monthly PPTA variance at the regional scale (r2 = 0.092). Although this proportion is modest in absolute terms, it constitutes a robust and physically interpretable signal within the highly variable tropical monsoon system, where local mesoscale convection, soil moisture feedback, and upstream Mekong discharge collectively influence precipitation variability [15,16]. Notably, the variance explained by ENSO increases substantially during the dry season transition months (November to April), with the Niño 3.4–PPTA correlation at lag 0 reaching −0.323 (p < 0.001), compared to −0.149 during the wet season (p < 0.05). This seasonal asymmetry demonstrates that ENSO impacts are most pronounced during monsoon withdrawal and the dry season in Southeast Asia, when background atmospheric moisture convergence is weaker and thus more susceptible to modulation by Pacific SST anomalies [15].

4.2. IOD Teleconnection: A Lagged Secondary Signal

The Dipole Mode Index demonstrates a teleconnection structure that differs markedly from that of Niño 3.4. Specifically, it shows weak and statistically insignificant correlations at short lags (0–7 months), but significant positive correlations at lags of 11–12 months (r = 0.186–0.182, p < 0.001). This long-lag positive relationship, which accounts for approximately 3.5% of the PPTA variance, warrants careful mechanistic interpretation.

The positive correlation at lags of 11–12 months aligns with the Indian Ocean Dipole's characteristic seasonal locking. Positive IOD events typically develop during June to October and dissipate by December [3,8]. A positive Dipole Mode Index during boreal summer (June–October) indicates anomalously warm sea surface temperatures in the western Indian Ocean and cool temperatures in the eastern Indian Ocean. This pattern suppresses moisture transport toward Southeast Asia through a weakened Somali Jet and anomalous anticyclonic circulation over the Bay of Bengal [10]. In the subsequent wet season (approximately 11–12 months later), residual thermocline adjustments in the eastern Indian Ocean may transition toward a negative IOD configuration, thereby enhancing eastward moisture flux and increasing VMD precipitation (Ashok et al., n.d.). This lagged recharge–discharge mechanism in the Indian Ocean thermocline is analogous to the delayed oscillator mechanism proposed for ENSO [35] and offers a physically plausible explanation for the observed 11–12 month lag.

The generally weak relationship between the Dipole Mode Index (DMI) and precipitation anomalies (PPTA) at the regional scale may be partially attributed to cancellation effects. Although positive Indian Ocean Dipole (IOD) events are typically associated with reduced VMD precipitation, these events often coincide with El Niño occurrences [10,13], complicating the separation of their individual impacts using a simple linear correlation. In the dataset, compound El Niño and positive IOD years, specifically 1997 (PPTA = +28.8%, Niño 3.4 = +1.11, DMI = +0.42) and 2023 (PPTA = +33.1%, Niño 3.4 = +0.82, DMI = +0.60), both displayed positive (above-normal) annual PPTA anomalies. This pattern suggests that in these instances, the subsequent La Niña transition (1998 and 2023–2024) may have exerted a dominant influence. Isolating the specific contribution of the IOD, independent of the ENSO-correlated component, would require partial wavelet coherence analysis or multivariate regression, which is proposed as a direction for future studies.

The spatial gradient of the DMI sensitivity was pronounced. The long-lag DMI signal was strongest at the upper delta and eastern coastal stations (r = 0.192–0.195 at lag 11) and weakest at Ca Mau (r = 0.140 at lag 12). These findings indicate that Indian Ocean teleconnections primarily affect the VMD through the South China Sea and Bay of Bengal moisture pathways, which exert a greater influence on the northeastern delta margins than on the western Gulf of Thailand coast [9].

4.3. PDO as a Decadal Background Modulator

The PDO demonstrated the most persistent and temporally stable teleconnection with the VMD PPTA among the three indices, as indicated by statistically significant negative correlations across all 12 lag values (r = −0.094 to −0.253, all p < 0.05). The PDO accounts for 6.4% of PPTA variance at the optimal lag of two months, a proportion comparable to the ENSO contribution (9.2%) and notably greater than the IOD contribution (3.5%). The nature of the PDO's influence is distinct: instead of acting through a specific atmospheric teleconnection pathway at a defined lag, it functions as a low-frequency background state that modulates the frequency and amplitude of ENSO events [11,12,14].

The phase-stratified annual PPTA analysis supports this interpretation in the following way. Years characterized by a positive PDO (PDO > 0.5) had a mean annual PPTA of −4.6%, whereas years with a negative PDO (PDO < −0.5) had a mean of +14.5%, representing a difference of approximately 19 percentage points. The observed transition in the VMD's precipitation regime around 2016, shifting from a persistent below-normal phase (2001–2015, predominantly positive PDO) to an above-normal phase (2016–2025, predominantly negative PDO), aligns with the documented 1998–2002 PDO phase shift toward negative values [12], highlighting the critical role of the PDO phase in establishing the decadal baseline for VMD water availability.

The modulation of El Niño impacts by the Pacific Decadal Oscillation (PDO) is clearly demonstrated in the historical record. El Niño events during positive PDO phases, such as those in 1982 (PPTA = −73.0%), 1987 (PPTA = −64.1%), and 2015 (PPTA = −26.7%), were linked to substantially more severe precipitation deficits than El Niño events during negative or neutral PDO phases, such as those in 2002 (PPTA = −18.7%) and 2023 (PPTA = +33.1%). The 1982–1983 El Niño, which coincided with a strong positive PDO (mean annual value of +1.17), resulted in the most extreme single-year deficit in the 45-year record (PPTA = −73.0%). This finding reinforces the concept of PDO-amplified ENSO forcing, as identified by [14] and [15] for the lower Mekong Basin.

The wavelet transform coherence (WTC) analysis supports this interpretation. Significant PDO–PPTA coherence is concentrated at 8–16 year periods, which are well above the 2–7 year ENSO band. The presence of significant coherence at both ENSO (2–5 year) and PDO (8–16 year) periodicities in the VMD precipitation signal indicates a hierarchical climate forcing structure. In this structure, the PDO establishes the decadal context within which ENSO-scale precipitation anomalies occur [13]. These findings have important practical implications. Seasonal climate forecasts for the VMD that consider only the ENSO state may be less effective during PDO phase transitions. Additionally, medium-range (5–10 year) water resource planning should explicitly incorporate the current PDO phase as a conditional baseline.

4.4. Spatial Gradient Across Hydrological Clusters

A consistent spatial gradient in teleconnection strength was evident across the three hydrological clusters. The southern–western coastal cluster (mean Niño 3.4 r = −0.287 at lag 2–3) and the eastern coastal cluster (mean r = −0.290) demonstrated slightly greater ENSO sensitivity than the upper delta cluster (mean r = −0.275). Although these differences are modest in magnitude, they remain physically meaningful and have distinct implications for water management throughout the delta.

The slightly stronger ENSO response observed at coastal stations, particularly at Ca Mau (r = −0.302), likely reflects the greater maritime influence at these sites. Both the East Sea and the Gulf of Thailand act as proximal moisture sources, and their sea surface temperature conditions are more directly coupled to Pacific ENSO forcing than the continental moisture transported via the Mekong River system to the upper delta [15]. In contrast, the upper delta stations (Chau Doc, Cao Lanh, and Can Tho) receive a substantial proportion of wet-season precipitation from the Mekong River flood pulse and associated moisture recycling, which partially decouples local precipitation from large-scale sea surface temperature forcing [16,21].

The 1-month difference in the optimal lag between the upper delta (lag 2) and southern coastal (lag 3) stations for Niño 3.4 indicates that ENSO signals reach the southern coastal periphery more slowly, likely due to Gulf of Thailand SST adjustments, compared to the direct atmospheric pathway to the upper delta. For the DMI, the 1-month lag difference between the eastern (lag 11) and southern coastal (lag 12) stations reflects varying exposure to the two main Indian Ocean moisture corridors: the Bay of Bengal pathway, which affects the northeastern delta via the South China Sea, and the Arabian Sea pathway, which influences the southwestern coast via the Gulf of Thailand [8].

These cluster-level differences have significant practical implications. Stations in the southwestern coastal cluster, which are already the most drought-vulnerable due to high dry-season saltwater intrusion and reliance on local rainfall for freshwater recharge, exhibited the strongest sensitivity to ENSO and IOD. Consequently, simultaneous El Niño and positive IOD events present the highest combined risk to freshwater availability and aquaculture productivity in the Ca Mau, Bac Lieu, and Rach Gia coastal provinces of Vietnam. In contrast, the upper delta is relatively buffered by upstream Mekong discharge; therefore, ENSO-induced precipitation deficits may be partially offset by sustained river inflows from the wetter upper Mekong headwaters in Yunnan and Tibet, where ENSO impacts differ from those in the delta [15,16].

4.5. Non-Stationarity and Compound Climate Forcing

The WTC analysis demonstrates the non-stationary nature of all three teleconnections. Significant coherence regions are episodic rather than continuous, and their frequency-band characteristics change throughout the 45-year record. This observed non-stationarity aligns with previously documented interdecadal modulations of the ENSO–monsoon teleconnection across Asia [13,30], highlighting the limitations of time-averaged Pearson correlation as a comprehensive measure of climate–precipitation linkages. The WTC results indicate that ENSO–VMD coherence at 2–5 year periods was especially pronounced during 1985–2000, which coincided with the active PDO+ phase and several major ENSO events, and again during 2014–2020. In contrast, coherence weakened substantially during 2000–2012, a period characterized by a more neutral PDO phase and relatively moderate ENSO events.

Compound event analysis provides additional evidence of non-linear interactions among climate modes. For example, the 1988 La Niña event, which coincided with a negative Indian Ocean Dipole (IOD; DMI = −0.201), resulted in the third-most extreme negative precipitation anomaly (PPTA = −64.3%) on record, despite the typical expectation of increased precipitation during La Niña events. This paradox likely reflects the strong negative IOD's suppression of Indian Ocean moisture flux toward Southeast Asia, which partially offsets the precipitation enhancement usually associated with La Niña [3,10]. In contrast, the 1996 La Niña combined with a negative IOD produced an anomalously wet year (PPTA = +42.8%), consistent with the additive precipitation-enhancing effects of the negative ENSO and negative IOD through both Pacific and Indian Ocean pathways. These findings underscore the limitations of analyzing teleconnection indices in isolation from other factors. Therefore, future operational drought and flood forecasting for the VMD should adopt a compound climate index framework that simultaneously monitors the ENSO phase, IOD state, and PDO background rather than relying on a single index. Developing such a compound index, weighted by the lag-specific correlation strengths identified in this study, offers a promising approach for improving seasonal precipitation outlooks for the delta [10,36].

Table 5 places the findings of this study in the context of key prior investigations on ENSO and Indian Ocean teleconnections with Southeast Asian precipitation. The results are broadly consistent with the regional literature but extend it in three important ways.

4.6. Implications for Water Resource Management and Climate Adaptation

The teleconnection signals identified in this study provide actionable insights for climate adaptation and water resource management in the VMD. The 2–3 month ENSO lag offers a practical window for seasonal drought and flood early warnings. Monitoring Niño 3.4 conditions in March and April enables a probabilistic forecast of the May–November wet-season precipitation anomaly, which is the primary determinant of annual rice production, aquaculture productivity, and reservoir storage throughout the delta.

The PDO phase should be incorporated as a conditional state variable in decadal water resource planning. The current negative PDO phase, which has persisted since approximately 2014 to 2016, is associated with a 14.5% increase in the mean annual PPTA anomaly in the dataset. This pattern indicates that VMD water availability may remain above the long-term average in the near future. However, PDO phase transitions, which have historically occurred every 20 to 30 years, could shift the baseline toward drier conditions, as observed during the positive PDO phase from 1977 to 1998. Infrastructure investment decisions, including reservoir capacity, irrigation channel design, and flood retention areas, should be stress-tested against both PDO phase scenarios.

Third, compound event risk, particularly the co-occurrence of El Niño and a positive IOD, constitutes the highest-impact scenario for the southwestern coastal cluster. The 1982–1983 (PPTA = −73.0%), 1987 (PPTA = −64.1%), and 2015–2016 compound events resulted in multi-season precipitation deficits that severely stressed agricultural water supplies and intensified saline intrusion in the Ca Mau Peninsula and Bac Lieu and Rach Gia coastal provinces. Concurrent monitoring of Niño 3.4 and DMI values, together with the lag structures identified in this study, should be incorporated into an operational drought risk assessment framework.

4.7. Limitations and Future Research Directions

This study has several limitations. Although the station-based precipitation dataset spanned a continuous 45-year period, it included only ten stations distributed across the VMD. This spatial density is likely insufficient to capture fine-scale precipitation gradients, particularly in the complex estuarine zones of the eastern delta, where tidal forcing and sea breeze convection generate highly localized precipitation patterns. Future research utilizing high-resolution gridded datasets (such as CHIRPS or PERSIANN-CDR) or denser rain gauge networks could enable more rigorous spatial interpolation and improved cluster delineations.

The analysis was limited to linear Pearson correlation and wavelet transform coherence, which do not capture non-linear or threshold-dependent relationships between climate indices and VMD precipitation. The severity-stratified analysis demonstrates that strong ENSO events produce non-linear responses that are inadequately represented by linear correlation. Employing non-linear methods, such as generalized additive models, quantile regression, or machine learning approaches, could provide a more comprehensive characterization of the ENSO–precipitation relationship across the entire range of index values.

While this study identifies the contributions of Niño 3.4, DMI, and PDO individually, the multi-variate interactions between these indices, including their partial (conditional) contributions after controlling for ENSO–IOD co-occurrence, were not formally quantified. Partial wavelet coherence analysis or multiple linear regression with variance inflation factor assessment would allow the independent contributions of each mode to be rigorously disentangled from their correlated components.

The analysis does not account for other potentially relevant climate forcing agents, including the Madden-Julian Oscillation (MJO), Boreal Summer Intraseasonal Oscillation (BSISO), Western North Pacific Subtropical High (WNPSH), or South China Sea SST anomalies, all of which can modulate seasonal precipitation in the VMD on sub-seasonal to interannual timescales. Incorporating these additional predictors into a multi-variate teleconnection framework represents an important priority for future work aimed at improving operational seasonal rainfall forecasting for the delta.

5. Conclusions

The spatiotemporal variability of monthly precipitation and its teleconnection with three large-scale climate indices (Niño 3.4, DMI, and PDO) were analyzed across ten meteorological stations in the VMD from 1981 to 2025 using Pearson lag-correlation analysis and the WTC. The findings indicate that the El Niño–Southern Oscillation (ENSO) is the primary driver of interannual precipitation variability in the VMD, exhibiting a peak negative correlation at a 2-month lag (r = −0.304, p < 0.001), accounting for approximately 9.2% of the regional monthly precipitation anomaly (PPTA) variance. This 2–3 month lag is shorter than those previously reported for inland Southeast Asian locations. The ENSO teleconnection was most pronounced in the southern and western coastal clusters, particularly during the dry season transition months, when background atmospheric moisture convergence was at its weakest.

The IOD exerted a secondary yet statistically significant lagged influence on VMD precipitation, with the strongest correlation observed at a lag of 11 to 12 months (r = 0.186, p < 0.001). This finding aligns with the thermocline recharge and discharge mechanisms in the Indian Ocean that operate on an annual cycle. Although the Pacific Decadal Oscillation (PDO) is not a direct forcing agent, it functions as a persistent low-frequency background modulator. Positive PDO phases reduced regional annual precipitation by an average of 4.6% relative to the long-term mean, whereas negative PDO phases increased it by 14.5%. The significant wavelet transform coherence was concentrated at periods of 8 to 16 years. Collectively, the three indices accounted for approximately 19% of the total monthly precipitation anomaly (PPTA) variance, with contributions from ENSO, PDO, and the DMI at 9.2%, 6.4%, and 3.5%, respectively.

The WTC analysis demonstrated that all three teleconnections exhibited non-stationarity, with significant coherence regions occurring episodically and concentrated within specific frequency bands. The most severe annual precipitation deficits in the record, such as those in 1982 (PPTA = −73.0%) and 1987 (PPTA = −64.1%), coincided with compound El Niño and positive PDO conditions. In contrast, simultaneous La Niña and negative IOD events were associated with the largest positive precipitation anomalies. The severity-stratified Pearson correlation analysis revealed a non-linear relationship between ENSO and PPTA, with the strongest correlations observed during weak-to-moderate ENSO phases rather than during extreme events. This finding underscores the limitations of linear indices in forecasting precipitation anomalies during major El Niño or La Niña years.

These findings have significant implications for seasonal climate forecasting and water resource management in the VMD region. The 2–3 month lead time associated with ENSO offers a practical prediction window for wet-season precipitation outlooks, whereas the current negative PDO phase indicates near-term projections of above-normal water availability. Nevertheless, the combined risk from simultaneous El Niño and positive IOD events, especially for drought-prone and saline-intrusion-vulnerable southern and western coastal provinces, underscores the need to integrate multi-index monitoring into operational early warning systems. Future research should employ partial wavelet coherence analysis to clarify the independent contributions of each climate mode, incorporate additional forcing agents such as the Madden-Julian oscillation, and expand the analysis to gridded satellite precipitation datasets to ensure spatially continuous coverage across the delta.