Marine heatwaves (MHWs), defined as anomalously warm sea surface temperature (SST) events, are analyzed over the Indian Ocean using NOAA OISST V2 (0.25° resolution) from 1981 to 2022, focusing on frequency, heatdays, maximum intensity, and mean intensity. Spatial analysis reveals MHW hotspots in the central-southwestern Indian Ocean, Bay of Bengal, and northern Arabian Sea, with 1–3 events per year lasting 3–15 days. Basin-wide mean maximum and mean intensities are 1.1 °C and 0.9 °C, respectively. Significant upward trends are found in all MHW metrics, including frequency (0.73 events/decade), heatdays (1.7 days/decade), mean intensity (0.11 °C/decade), and maximum intensity (0.14 °C/decade), with the Arabian Sea exhibiting the largest rise in heatdays (17.79 days/decade). A comparison of the periods 1981–2000 and 2001–2022 shows increased frequency and duration in the latter, though with marginally lower intensities. Empirical Orthogonal Function (EOF) analysis explains 28% and 17% of the spatial variance in frequency and duration, respectively. Partial correlation and kernel-based methods identify the El Niño–Southern Oscillation (ENSO), particularly positive Niño3.4 anomalies, as the dominant driver of MHW variability, strongly correlated with frequency (r = 0.59) and intensity (r = 0.5). Granger causality and Convergent Cross Mapping confirm ENSO’s particularly Niño3.4’s causal influence and reveal bidirectional feedback with MHWs heatdays, linked to Walker circulation driven changes in zonal winds and thermocline depth. The Pacific–Japan (PJ) pattern exerts secondary influence, while the Indian Ocean Dipole (IOD) plays a limited role. This underscores the importance of incorporating such drivers in predictive models and regional adaptation strategies.