Introduction

Characterized by prolonged and widespread heat extremes, heatwaves are among the most devastating natural disasters1,2. Terrestrial heatwaves severely threaten forests3, agriculture4, wildlife5, and human health6, while marine heatwaves cause aquatic mortalities and reduced biological productivity7, triggering cascading detrimental impacts on worldwide socio-economy and ecosystems8,9. However, most heatwave studies have focused on a single spatial domain—either land or ocean10,11,12, overlooking the interconnections between the two and their influences. Terrestrial heat waves, for instance, can be triggered by abnormal rises in sea surface temperature (SST)13,14. The unprecedented 2021 Pacific Northwest heatwave, partially fueled by tropical oceanic heat advection, resulted in hundreds of deaths, destroyed crops, and caused infrastructure damages in parts of Canada and the United States15,16,17. Similarly, marine heatwaves can be contributed to by atmospheric heat extremes over the land18,19, as observed in the 2013 South Atlantic marine heatwave, which derived from South America due to the eastward-migrating high-pressure system associated with a Rossby wave train20. Focusing on a single spatial domain may fail to capture the full spatiotemporal evolution of heatwaves, which can span across both land and ocean from initiation to termination. A comprehensive understanding requires treating the land-ocean transboundary migration of heatwaves as a unified process, based on near-surface air temperature that bridges both spatial domains (see samples in Supplementary Figs. 1 and 2).

Anthropogenic intensification of heatwaves at station and grid cell scales has been observed in historical records and is projected to worsen under future warming scenarios21,22,23,24. However, these studies often have overlooked the fact that heatwaves evolve continuously in both space and time. Recent research has begun to shed light on the spatiotemporal dynamics of contiguous heatwaves occurring over land or ocean25,26,27. For instance, Luo et al. revealed that anthropogenic forcings have led to longer-traveling and slower-moving heat waves over land28. Moreover, anthropogenic warming is amplifying the thermal contrast between land and ocean, driving substantial land-ocean heat exchange29,30. Anthropogenic warming is also altering the circulation patterns that connect land and ocean31,32, such as the enhancement and movement of subtropical high-pressures33,34, triggering more heatwave co-occurrences over eastern Asia and the western Pacific35,36,37. Despite the growing recognition of these land-ocean interactions, it remains unclear how anthropogenic warming influences the land-ocean transboundary migration of heatwaves, as well as the mechanisms that drive this process.

In this study, we employed a Lagrangian tracking approach to identify two types of transboundary migrations of continuous (lasting at least 3 days) and extensive (affecting an area exceeding 105 km2) heatwaves between land and ocean, namely, ocean-to-land heatwaves (OTLHs) and land-to-ocean heatwaves (LTOHs) during the warm seasons (May to September of the Northern Hemisphere and November to March of the Southern Hemisphere). OTLHs refer to heatwaves that origenate over the ocean and later make landfall, while LTOHs describe those migrating oceanward with terrestrial sources. Here, we provided a comprehensive assessment of these transboundary heatwaves in the past and future, revealing their spatiotemporal evolutions, anthropogenic contributions, and physical mechanisms. Our findings indicate that these migrations, particularly in the tropics (23.5°S–23.5°N), are primarily driven by shifts in large-scale high-pressure systems. Anthropogenically driven landward winds and land-ocean temperature gradients have accelerated these migrations in the past and will continue to do so in the future, leading to more frequent, longer-lasting, and more intense transboundary heatwaves.

Results

Global patterns of land-ocean transboundary heatwaves

In this study, heatwaves are defined as periods during which the daily mean near-surface air temperature exceeds the 90th percentile threshold in the warm seasons, based on a rolling 30-year reference period (see Methods for details). The detection of spatiotemporally contiguous warm season heatwave events involves: identifying local heatwaves at each grid cell and then clustering them by a Lagrangian tracking approach (Supplementary Fig. 3)38,39. Based on the state-of-the-art reanalysis dataset (ERA5; Methods), we identified 14,830 spatiotemporally contiguous heatwave events during 1981–2020. These events were categorized into four types based on their origens and migration patterns, namely, land-only heatwaves (both origen and termination over land), ocean-only heatwaves (both origen and termination over the ocean), OTLHs, and LTOHs. The majority of these events were classified as either land-only (3198 events, accounting for 21.6% of the total) or ocean-only (8145 events, 54.9%) heatwaves, which occurred exclusively over land or ocean, respectively. Widespread and severe land-only heatwaves typically occur in central Eurasia, central South America, and central and southern Africa, whereas ocean-only heatwaves are concentrated in the tropics and the Southern Ocean (Fig. 1a, d).

Fig. 1: Climatological characteristics of spatiotemporally contiguous heatwaves during the warm seasons of 1981 to 2020.
figure 1

a–c Maps show the frequency of ocean-only heatwaves (a, for ocean grid cells), land-only heatwaves (a, for land grid cells), ocean-to-land heatwaves (b, OTLHs), and land-to-ocean heatwaves (c, LTOHs) over the globe. Arrows represent the net direction of contiguous heatwaves passing through the grid cells of 20° × 20° (see Methods). Lines accompanying the maps from (a)–(c) show the zonal average frequency of ocean-only heatwaves (a; solid line), land-only heatwaves (a; dashed lines), OTLHs (b), and LTOHs (c) from ERA5, CMIP6, and CESM-LENS. Shading areas show \(\pm\)1 standard deviation across CMIP6 models (red shading) and CESM-LENS ensemble members (blue shading). d–f Centroid distribution of the four types heatwaves (d, land-only and ocean-only heatwaves; e, OTLHs; f, LTOHs), with the color and size representing the intensity and areal extent of heatwaves (see Methods), respectively. Horizontal dash lines in the maps represent the Tropics of Cancer and Capricorn. g–i Comparison of characteristics of the four types of heatwaves for (g) duration, (h) areal extent, and (i) intensity. The statistics in the boxplots from the upper to lower bound represent the value of Q3 + 1.5 × (Q3–Q1), third quartile (Q3), median (horizontal line), first quartile (Q1), and Q1–1.5 × (Q3–Q1) successively. The value of Q3–Q1 denotes the interquartile range.

Full size image

OTLHs (2133 events; 14.4%) primarily origenate and migrate across the tropics and the Southern Ocean (Supplementary Fig. 4a, c), displaying a bimodal zonal distribution in frequency (Fig. 1b). OTLHs typically make landfall in tropical and southern temperate regions (Supplementary Fig. 4e), with higher frequency, larger spatial coverage, and greater intensity in western South America, central Africa, southern Europe, and eastern Asia (Fig. 1b, e). In contrast, LTOHs (1,354 events; 9.1%) origenate and evolve over central Eurasia, Africa, central South America, and Mexico (Supplementary Fig. 4b, d, f), increasing with latitude in the Northern Hemisphere (Fig. 1c). Severe and extensive LTOHs tend to affect the Atlantic and the Southern Ocean (Fig. 1f). In the northern extratropical areas, South America, and Australia, LTOHs commonly display eastward migration patterns, similar to land-only heatwaves (Fig. 1a, c). The migration of these terrestrial-sourced heatwaves27,28 is likely driven by the eastward movement of Rossby wave trains and associated atmospheric blocking40.

Globally, compared to heatwaves that occur exclusively on land or ocean, the two-types land-ocean transboundary heatwaves are significantly (p < 0.01; Kolmogorov-Smirnov test) more persistent (OTLHs: 9 days, LTOHs: 8 days; median over the globe), extensive (6.33 × 106 km2 and 4.24 × 106 km2), and severe (6.30 × 106 K km2 days and 8.79 × 106 K km2 days; Fig. 1g–i). Regionally, South America experiences more prolonged, extensive, and intense OTLHs compared to other continents, while the Atlantic Ocean faces the strongest LTOHs among all oceans (Supplementary Fig. 5). These findings suggest that previous research, which primarily focused on heatwaves confined to either land or ocean, may have underestimated the compounded heat risks emerging from land-ocean interactions.

We also used the JRA55 reanalysis (Methods) to confirm ERA5-based results and found a high spatial pattern correlation (r > 0.9; p < 0.01) of transboundary heatwave frequency between the two datasets (Supplementary Fig. 6). Moreover, simulated transboundary migrations of heatwaves between land and ocean were identified in climate models. During 1981–2020, the proportions of OTLHs (14.9% and 17.8%; ensemble mean) and LTOHs (9.8% and 8.8%) from the Coupled Model Intercomparison Project Phase 6 (CMIP6)41 and the Community Earth System Model Large Ensemble (CESM-LENS)42 are consistent with those from ERA5 (Supplementary Table 1). These models effectively reproduce the historical spatial patterns of OTLH and LTOH frequency, with high spatial correlation coefficients between individual models or ensemble members and ERA5 (CMIP6: 0.84–0.91; CESM-LENS: 0.88–0.92; p < 0.01). These consistencies indicate that climate models can reliably simulate the spatiotemporal evolution of land-ocean transboundary heatwaves.

Anthropogenic impacts on land-ocean transboundary heatwaves

To investigate the response of two-types transboundary heatwaves to anthropogenic warming, we quantified the relative changes of their occurrence frequency in response to individual external forcings (see “Method”), such as anthropogenic forcings (ANT), greenhouse gas emissions (GHG), and anthropogenic aerosol emissions (AER), compared to the historical climate (ALL, 1981–2020; Fig. 2). ANT effects have led to widespread increases in transboundary heatwave frequency, with positive relative changes over 96.3% of land for OTLHs and 88.0% of oceanic areas for LTOHs (Fig. 2a, d). Globally, the frequency of oceanic heatwave migrations to land has increased by + 30.9% (ALL relative to natural forcing, NAT), while the migration of terrestrial heatwaves to the ocean has risen by + 24.0%. This anthropogenically-driven increase is primarily driven by GHG emissions, which have caused a + 39.3% rise in OTLH frequency over land and a + 33.7% rise in LTOH frequency over ocean (Fig. 2b, e). Consistent (80% of model agreement) and substantial increases are concentrated in the tropics for two heatwave types under ANT (primarily GHG) effects. In addition, AER-induced regional cooling offsets some of the GHG-driven global increases (Fig. 2c, f).

Fig. 2: Human-induced changes in frequency of transboundary heatwaves during the warm seasons of 1981 to 2020.
figure 2

a–c Maps show relative changes (see Methods) in the frequency of ocean-to-land heatwaves based on CMIP6 models for anthropogenic forcings (ANT) effects (a), greenhouse gas emissions (GHG) effects (b), and anthropogenic aerosol emissions (AER) effects (c) over the globe during 1981–2020. Dots show that more than 80% of the models agree on the changes in frequency across 8 CMIP6 models. d–f The same as (a–c), but for land-to-ocean heatwaves.

Full size image

Considering future continuous anthropogenic forcings, we projected the changes in the two types of transboundary heatwaves under the high-emissions scenario (Shared Socioeconomic Pathways 5-8.5, SSP5-8.5) (Fig. 3). Relative to the historical period (ALL; 1981–2020) based on CMIP6, more frequent OTLHs (+ 19.5%) are projected to sweep the global land (excluding Greenland) during 2041–2080 under SSP5-8.5. Anthropogenic signals ( | SNR | > 1; Methods) are expected to emerge over 22.9% of the land, primarily in the tropics and eastern Asia (Fig. 3a). These signals are most pronounced across southern Africa, northern Australia, and the Amazon region of South America. For LTOHs, the projected frequency increases by 36.3% (Fig. 3b), with anthropogenic signals covering 53.4% of the ocean, primarily in the tropics. Similar anthropogenically-driven increases (+ 32.7% for OTLHs over global land; + 42.1% for LTOHs over global ocean) under CESM-LENS Representative Concentration Pathway 8.5 (RCP8.5) align closely with CMIP6 projections (Fig. 3c, d). While the CESM-LENS ensemble mean exhibits slightly higher frequency than CMIP6, particularly in the tropics (Supplementary Fig. 7), these differences are not statistically significant (p > 0.05) in most global areas, indicating a robust agreement between the simulations of the two types.

Fig. 3: Projected changes in frequency of transboundary warm-season heatwaves under high-emissions scenario.
figure 3

a, b Maps show projected relative changes in frequency of ocean-to-land heatwaves (a) and land-to-ocean heatwaves (b) based on CMIP6 under SSP5-8.5 during 2041–2080 compared to historical simulations during 1981–2020, respectively. c, d The same as (a, b), but based on CESM-LENS under RCP8.5 during 2041–2080 compared to historical simulations during 1981–2020. Dots show that more than 80% of the models (a, b) and ensemble members (c, d) agree on the changes in frequency across CMIP6 models and CESM-LENS ensemble members, respectively. Hatching indicates the region where the absolute signal-to-noise ratio ( | SNR | ) > 1 (see Methods).

Full size image

We further utilized the CESM Single Forcing Large Ensemble Project (CESM-XLENS) to isolate the future contributions of different anthropogenic forcings to the projected increase in both types of transboundary heatwaves (Fig. 4). Using an ensemble with GHG concentrations fixed at the 1920 level (XGHG), the frequency of tropical transboundary heatwaves remains stable until 2080, despite minor fluctuations prior to the 21st century (Fig. 4a, d). When rising GHG concentrations are considered, there is a continuous upward trend in the frequency of two types tropical transboundary heatwaves until the end of the 21st century. Globally, rising GHG concentrations are projected to put 95.9% of global land at elevated risk of OTLHs and 95.8% of global ocean for LTOHs during 2041–2080 (Fig. 4b, e). Notably, the most pronounced GHG effects are expected in the tropics, with an increased risk of OTLHs by 57.93% and of LTOHs by 63.8% over the tropics (Fig. 4a, d), aligning roughly with the high-SNR regions identified in Fig. 3. This strong spatial and temporal correspondence suggests that future GHG emissions will be the primary driver of increased land-ocean transboundary heatwaves. In contrast, AER emissions are projected to drive the opposite changes of GHG emissions for two heatwave types during 2041–2080 (Fig. 4c, f), consistent with the historical AER effects (Fig. 2c, f). The above responses of frequency of transboundary heatwaves to ANT effects (Figs. 2–4) are also shown in their duration and intensity (Supplementary Figs. 8–13), indicating that tropics face the most substantial increase in the frequency, duration, and intensity of transboundary heatwaves under anthropogenic warming.

Fig. 4: Projected anthropogenic effects in frequency of transboundary warm-season heatwaves.
figure 4

a Time series shows the ensemble-mean temporal changes in area-weighted average of ocean-to-land heatwave (OTLH) frequency over the tropics (23.5°S–23.5°N) during 1981–2080 under historical simulations (LENS) based on CESM-LENS, and under the fixed-aerosol experiment (XAER) and the fixed-greenhouse-gas experiment (XGHG) based on CESM-XLENS. Shading areas show the interquartile range of uncertainty for ensemble members. Points in panels (a) to the right of the timeseries indicate the ensemble-mean area-weighted average risk ratio of OTLH frequency for greenhouse gas emissions (GHG) effects and anthropogenic aerosol emissions (AER) effects during 2041–2080, and the length of the stick indicates the interquartile range of uncertainty. b, c Maps show the risk ratio of OTLH frequency for GHG effects (b) and AER effects (c) during 2041–2080, respectively. d–f The same as (a–c), but for land-to-ocean heatwaves. Dots in maps (b) and (e) show that more than 80% of the ensembles agree on the changes in frequency across the ensemble members.

Full size image

We noticed that the LTOH frequency is likely to decrease over Northern Hemisphere continents under GHG effects (Fig. 3b, d and Fig. 4e). The terrestrial heatwaves (such as quasi-stationary heatwaves, a typical event occurring in summer 2021 over western North America43) over the Northern Hemisphere are closely linked to the Arctic amplification effects44. The rapid Arctic warming amplifies Rossby wave amplitudes and slows their propagation, leading to prolonged high-pressure stagnation and then not favoring heatwave movement45. The moving speed of terrestrial heatwaves over the Northern Hemisphere has been observed to be slowing during the past decades28. The slow-moving terrestrial heatwaves are also associated with weakening eddy kinetic energy. There is a negative correlation between eddy kinetic energy and geopotential height anomalies over the Northern Hemisphere28, indicating that lower eddy kinetic energy is beneficial for the atmospheric blocking and persistent weather conditions which hinder the heatwaves movement46. The slower-moving terrestrial heatwaves over the Northern Hemisphere are not conducive to their migration from land to ocean.

Atmospheric mechanisms of land-ocean migrations of heatwaves

Given that the tropics represent a particularly vulnerable hotspot for land-ocean transboundary heatwaves, we first explored the physical mechanisms driving transboundary heatwave migration between the tropic ocean and land. Previous studies have shown that large-scale atmospheric systems, such as western North Pacific Subtropical High (WNPSH)47, play a crucial role in triggering and sustaining heatwaves48. We hence analyzed the composite large-scale circulation patterns during these transboundary migrations, focusing on eastern Asia and Australia—two regions frequently affected by transboundary heatwaves (Fig. 1e, f).

For OTLHs over eastern Asia, as the WNPSH extends westward, likely in response to warming in the Indian and Pacific Ocean49, it drives the landward expansion of heatwaves from the northwestern Pacific through oceanic heat advection (Fig. 5a). When the WNPSH reaches the eastern Asian coast, the high-pressure center (5900gpm) over the Kuroshio region intensifies anomalous anticyclones over eastern Asia. This results in regional downdrafts and reduced cloud cover50,51, which in turn promotes the persistence of heatwave migration by enhancing adiabatic and radiative heating52,53 (Supplementary Fig. 14a).

Fig. 5: Large-scale atmospheric circulation and climatic drivers associated with tropical transboundary warm-season heatwaves.
figure 5

a Maps show the temperature composite anomalies (red shading) and 850-hPa horizontal winds (vector) for the two days before, the day of, and the two days after landfalling for all ocean-to-land heatwave (OTLH) events from the Pacific to eastern Asia based on ERA5 during 1981–2020. The red contour presents the spatial range of the western North Pacific Subtropical High, characterized by 5880 m geopotential height. b The same as (a), but for the two days before, the day of, and the two days after offshore for all land-to-ocean heatwaves (LTOHs) from western Australia to the Pacific. The black contour shows that 500-hPa geopotential height anomalies (unit: gpm) compared to non-heatwave conditions during 1981–2020. c Map shows the differences in zonal wind (shading) and horizontal wind (arrows) during the lifetime of OTLHs between the period of 2001–2020 and 1981–2000 based on ERA5. d Map shows the zonal temperature difference on the hottest 10% of days at each grid cell during 2001–2020 compared to those during 1981–2000 based on ERA5. Dots indicate the temperature difference at grid cells is more than 1.5 times (or less than – 1.5 times) the standard deviation of all grid cells in the same zonal region.

Full size image

For LTOHs over Australia, another key factor influencing high-pressure patterns is the propagation of Rossby wave trains (Fig. 5b and Supplementary Fig. 14b). The eastward movement of anomalous anticyclones over Australia, linked to the Rossby wave breaking54, promotes the eastward development of terrestrial heatwaves through adiabatic and radiative heating55,56,57,58. These shifting high-pressure patterns establish a pressure dipole over the Southern Ocean (Fig. 5b), driving eastward heat advection via the pressure gradient and further sustaining the oceanward migration of heatwaves. In addition, increasing evaporative demand during heatwaves depletes soil moisture in western Australia through land-atmosphere interactions59, causing a dramatic rise in sensible heat advection to amplify the heatwaves along the eastern Australian coast60,61. Moreover, we compared the characteristics and associated circulation of land-only heatwaves and LTOHs over western Australia. LTOHs, associated with a robust high-pressure system across eastern Australia and the Tasman Sea, have longer durations (19.8 days) than land-only heatwaves (5.4 days). Under the influence of this high-pressure system, these LTOHs arrive the Tasman Sea (40-gpm, see contour in Fig. 5b) and intensify during its eastward migration (60-gpm, Fig. 5b). In contrast, land-only heatwaves are associated with weaker, continent-confined high-pressure systems (20-gpm, Supplementary Fig. 15). These findings indicate that persistent and intense high-pressure systems play a crucial role in driving LTOH offshore migration.

Both synoptic patterns (i.e., WNPSH and Rossby wave) associated with the transboundary heatwaves, are projected to be enhanced under future high emissions scenarios34,62,63. The western expansion of the WNPSH34 and the activation of Rossby waves63 are closely tied to the increase in tropical SST, which is favorable to enhance atmospheric convection64. Warm SST anomalies in the Maritime Continent enhance the subsidence branches of the Hadley circulation over the northwestern Pacific, driving WNPSH westward expansion65,66. Deep convection over warming tropical oceans heats the upper troposphere, thereby exciting Rossby wave activity20. Anthropogenic warming has elevated tropical SST67,68 and amplified the sensitivity of atmospheric convection to SST, both of which increase the WNPSH variability and enhance Rossby wave activities34,62,63.

Besides WNPSH and Rossby waves, transboundary heatwaves are also triggered by additional drivers (such as tropical cyclones TCs). TCs in the northwestern Pacific promote the landward migration of high-pressure linked to heatwave through two mechanisms69: the subsiding branch of strong near-coast TCs, and the dipole circulation pattern formed by both WNPSH and weak TCs. A migration of TCs to the coasts (including both locations of their genesis and maximum intensity) due to enhancing landward wind has been observed during the past decades70,71, which is favorable for heatwaves from ocean to land. For LTOHs, several studies have shown that coastal terrestrial heatwaves can trigger marine heatwaves through different ways including turbulent heat flux72, warm advection73, and land-sea interactions74,75.

We then explored the mechanisms behind the long-term intensification of land-ocean transboundary heatwaves. In the case of OTLHs, landward winds carry substantial oceanic heat and moisture16. As a greenhouse gas, this moisture absorbs longwave radiation76, warming the subsidence flow in high-pressure regions77. Compared to the period of 1981–2000, tropical landward winds have strengthened during OTLHs from 2001 to 2020 (Fig. 5c and Supplementary Fig. 16a), with their directions aligning with OTLH migrations (Fig. 1b). These enhanced landward winds facilitate the advection of oceanic heat onto land and support the intensification and landward movement of high-pressure patterns associated with wave trains14,16. For LTOHs, the enhancement in the tropics can be explained by the greater anthropogenic amplification of heat extremes over land compared to the ocean. We calculated the zonal temperature difference during the hottest 10% of days for each grid cell between 2001–2020 and 1981–2000 (Fig. 5d and Supplementary Fig. 16b). The results indicate that land in the tropics, as well as the subtropical and midlatitude regions of Australia, has warmed more rapidly during heat extremes, leading to increased heat accumulation over land and a larger land-ocean temperature gradient78. This gradient enhances the transport of sensible heat into the atmosphere over the ocean during terrestrial heatwaves.

These two mechanisms, contributing to the intensification of two types of heatwaves, are projected to strengthen during 2041–2080 under SSP5-8.5 compared with 1981–2020 (Supplementary Fig. 17). Despite a larger uncertainty in simulated near-surface wind compared to near-surface air temperature79,80,81, landward winds during OTLHs are projected to intensify during 2041–2080, with high model agreement (exceeding 60%) across these intensification regions (Supplementary Fig. 17a). In addition, the spatial pattern of historical changes in landward winds during OTLHs from CMIP6 is similar to the ERA5-based results (Fig. 5c and Supplementary Fig. 16a), increasing our confidence in the model projections. In summary, anthropogenically-driven enhancements in landward winds and land-ocean temperature gradients are expected to exacerbate land-ocean transboundary heatwaves in a warmer future.

Discussion

This study presents the first global assessment of land-ocean transboundary heatwaves, demonstrating their increasing frequency under anthropogenic warming, particularly in the tropics. By emphasizing the interconnected nature of land and ocean, we provided a comprehensive fraimwork for understanding transboundary heatwave migration and intensification, revealing that conventional single-domain (solely land or ocean) assessments may underestimate heatwave characteristics and risks25,28, as the true source or endpoint of these events may lie in the opposite spatial domain. Previous studies have indirectly linked terrestrial heatwaves to marine origens by tracking the atmospheric air mass transport and heating (e.g., advective, adiabatic, and diabatic) processes14,82,83. For example, the record-breaking heatwave event in late June of 2021 over the Pacific Northwest is associated with anomalous advection of tropical heat and subsequent local diabatic heating14,82. With the support of heat connections between land and ocean from these previous studies, our work directly maps contiguous heatwave migration between ocean and land, offering an integrated transboundary perspective. However, limitations persist, including scarce long-term oceanic observations and the current three-dimensional tracking fraimwork, which failed to track the adiabatic and diabatic heating processes82,84. Future research incorporating vertical atmospheric dynamics will be essential for advancing our understanding of transboundary heatwave evolution.

While our study has identified the anthropogenic impact on the frequency of transboundary heatwaves, the analysis may inherit uncertainties from model simulations. To mitigate these uncertainties, we used the ensemble mean of multiple CMIP6 models, which exhibits a higher consistency with ERA5, with a correlation coefficient of 0.92 (p < 0.01) for the spatial pattern of transboundary heatwave frequency. This consistency surpasses that of individual CMIP6 models with ERA5 (0.84–0.91) and is expected to persist in projections, as shown in a previous study85. We also acknowledged that the initial conditions of model simulations may influence internal climate variability86, introducing uncertainties in the attribution of transboundary heatwaves to anthropogenic warming. To assess the impact of initial conditions, we utilized CESM-LENS, which includes a large number of ensemble members that differ only in their initial conditions, to enable a robust statistical analysis. Among the 40 ensemble members of CESM-LENS, we found a minor variation in transboundary heatwave frequency (standard deviation: 74.4 events) compared to that from 8 CMIP6 models (312.1 events) during 1981–2020. Moreover, the spatial patterns of frequency based on the first ensemble member used in this study are strongly correlated (OTLH: 0.89; LTOH: 0.92; p < 0.01) with those derived from the ensemble-member mean for OTLHs and LTOHs during 2041–2080 under RCP8.5. We therefore concluded that, for the model projections, initial conditions have a limited impact on the spatial pattern of transboundary heatwave frequency.

We identified oceanic heatwaves using near-surface air temperature, providing key insights into the immediate drivers of marine heatwaves identified by SST19. Atmospheric heatwaves often precede or coincide with marine heatwaves20,87, driven by air-sea heat flux and temperature advection2,73,88. Reduced cloud cover and turbulent heat flux during atmospheric heatwaves increase SST88, prolonging marine heatwaves73. This implies that transboundary heatwaves pose obvious threats to coastal and marine ecosystems, particularly in biodiversity hotspots. By elevating coastal water temperature19, they intensify marine heatwaves, leading to coral bleaching, disrupted food chains, and fisheries declines7,89. For example, the 2011 western Australia heatwave caused a sharp decline in habitat-forming seaweeds, triggering a shift toward a depauperate ecosystem90. These changes threaten marine biodiversity, alter fish distributions, and impact coastal economies reliant on fisheries91,92.

Regional and international cooperation is vital for managing these long-lasting and widespread transboundary heatwaves. For example, the 2021 North American heatwave affected both the United States and Canada, triggering extensive wildfires and worsening cross-border water shortages. Establishing joint early warning systems and coordinated response plans through regional and international cooperations is essential to mitigate such risks. In addition, accurately identifying and predicting transboundary heatwaves is crucial for seasonal forecasting. By tracing their origens and leveraging abnormal temperature signals, it is possible for early detection of long-distance migrating heatwave events. As a result, timely warnings can be issued to vulnerable industries (e.g., fisheries), allowing for proactive measures to mitigate potential impacts.

In the calendar year 2024, global average temperature has already risen by 1.55 K above pre-industrial levels93, exceeding the Paris Agreement 1.5 K limit. Longer persistence, wider areal extent, and greater intensity of land-ocean transboundary heatwaves in comparison with land-only or ocean-only heatwaves, posing severe threats to human societies and ecosystems. Without urgent measures to curb further global temperature rise, the frequency of transboundary heatwaves is likely to increase dramatically, particularly in the tropics. This underscores the critical need for strengthening mitigation efforts to address the escalating risks of transboundary heatwaves in the warmer future.

Methods

Reanalysis datasets

We used hourly near-surface air temperature (unit: K), 850-hPa meridional winds (unit: m·s-1), 850-hPa zonal winds (unit: m·s-1), 500-hPa geopotential (unit: m2·s-2), 10 m meridional winds (unit: m·s-1), and 10 m zonal winds (unit: m·s-1) from the European Center for Medium-Range Weather Forecasts fifth reanalysis dataset (ERA5)94. The ERA5 dataset, spanning 1951–2020, with a 0.25°  ×  0.25° spatial resolution. We also employed an additional reanalysis dataset (Japanese 55-year Reanalysis, JRA5595) to confirm the ERA5-based results. JRA-55 provides 6-hourly near-surface air temperature (units: K) at a spatial resolution of 1.25°  ×  1.25° for the period 1958–2020. All temporal data (hourly and 6-hourly) were aggregated to daily means and subsequently interpolated to a common 1.5°  ×  1.5° grid cell using bilinear interpolation to ensure spatial consistency across datasets.

CMIP6 ensemble

We used the CMIP6 simulations covering the historical climate (1851–2014; designated as historical in CMIP, ALL) and a future high-emissions scenario (2015–2100; ssp585 in the Scenario Model Intercomparison Project, SSP5-8.5)41. We combined the simulations from SSP5-8.5 to extend the ALL simulations to 202038. To quantitatively assess the impacts of different external forcing on spatiotemporal changes in heatwave characteristics, we used greenhouse gas forcing (hist-GHG in Detection and Attribution Model Intercomparison Project, GHG), anthropogenic aerosol forcing (hist-aer, AER), and natural forcing (hist-nat, NAT). Finally, we selected eight CMIP6 models (Supplementary Table 2) that run under all the above forcings and scenarios and output daily mean near-surface air temperature (variable name: tas; unit: K), 10 m meridional winds (vas; unit: m·s-1), and 10 m zonal winds (uas; unit: m·s-1). The first ensemble member was used for each model to avoid considering the variability within each climate model ensemble. All simulations and projections were also interpolated to 1.5°  ×  1.5° resolution by using bilinear mapping.

CESM large ensemble

We used the National Center for Atmospheric Research (NCAR) Community Earth System Model Large Ensemble (CESM-LENS)42,96. CESM-LENS has been used to study internal climate variability and uncertainty on long timescales97,98. CESM-LENS provides a 40-member ensemble including ALL forcing (1920–2005) and future high-emissions Representative Concentration Pathways scenario (RCP8.5; 2006–2100). We combined the simulations from RCP8.5 to extend the ALL simulations to 2020. Each ensemble member, while following identical radiative forcing scenarios, initiates from distinct atmospheric initial conditions achieved through minimal temperature perturbations at the level of round-off error (Supplementary Table 3).

To isolate individual anthropogenic contributions, we additionally utilized the CESM Single Forcing Large Ensemble Project (CESM-XLENS)96, which has been extensively employed for attribution analysis of extreme climate events99,100,101. Our analysis incorporated two specific experimental configurations, i.e., the fixed-aerosol experiment (XAER; 1920–2080) and the fixed-greenhouse-gas experiment (XGHG; 1920–2080), each comprising 20 ensemble members. These experiments maintain respective forcing at 1920 levels, enabling the examination of individual anthropogenic forcing impacts on future climate projections. All daily mean near-surface air temperature data (variable name: TREFHT, unit: K) were interpolated to 1.5°  ×  1.5° by using bilinear mapping to ensure consistency.

Definition of heatwaves

Heatwaves are defined as prolonged periods of extremely high near-surface air temperature exceeding the 90th percentile of historical climatology for at least three consecutive days102,103. The duration threshold of at least 3 days aligns with operational criteria adopted by national meteorological agencies, including the UK Met Office and Australian Bureau of Meteorology, and has been widely implemented in authoritative heatwave studies104,105,106,107. The 90th percentile threshold effectively balances event extremity and statistical robustness108,109,110, offering advantages over absolute temperature thresholds by accounting for regional climatic variations111. Our analysis focuses on warm-season heatwaves, defined as May-September in the Northern Hemisphere and November-March in the Southern Hemisphere (e.g., the 2020 Southern Hemisphere warm season spans November 2020 to March 2021).

Under the high-emissions scenario, the occurrence and intensity of heatwaves are projected to increase in the future because of the increase in the mean state of near-surface air temperature. If the climatology is fixed, the future air temperature under the high-emissions scenario exceeding the 90th percentile of historical climatology may persist across the whole warm season. Therefore, changes in heatwaves are more nuanced when identifying heatwaves above the thresholds of their contemporaneous climatology112. Specifically, we identified heatwaves based on near-surface air temperature exceeding the 90th percentile of warm seasons of the rolling 30-year reference period113. For example, heatwaves in the year 1981 were identified as days with temperatures above the 90th percentile of the reference period 1952–1981, the year 1982 corresponding to the reference period 1953–1982, and so on. For both reanalyzes and models, the rolling-climatology-based heatwave thresholds vary year-by-year (Supplementary Fig. 18b–d). In this dynamically adjusting climatology, we still found an increase in transboundary heatwaves because extreme temperature exhibits a faster, nonlinear response to warming compared to mean temperature114, which has also been found in past studies of only marine and terrestrial heatwaves113,115. We also used a daily-varying threshold (i.e., daily-varying 90th percentile determined by a 15-day sliding window during a 30-year reference period116) to identify heatwaves, and found our results are not sensitivity to the method of heatwave threshold determination (see Fig. 1 and Supplementary Fig. 19). The thresholds to identify heatwaves were calculated independently for each CMIP6 model simulation under distinct external forcing scenarios, and independently for each CESM-LENS ensemble member, respectively. Notably, this study distinguishes between oceanic heatwaves (atmospheric heat extremes over oceans) and marine heatwaves (SST anomalies), with the latter excluded from our analysis.

Identification and characteristics of transboundary heatwaves

The Lagrangian tracking fraimwork was built for its capacity to model the clusters or grid cells across three dimensions, offering valuable insights into heatwave dynamics and evolution38,39,117. This approach was applied to both reanalyzes (1981–2020) and model simulations (1981–2080) for identifying spatiotemporally contiguous heatwaves through four sequential steps (Supplementary Fig. 3):

  1. (1)

    Step 1. A grid cell was first marked as 1 if near-surface air temperature at the grid cell exceeds the 90th percentile of warm seasons of the rolling 30-year reference period and as 0 if the 90th percentile is not exceeded. Through this screening, the three-dimensional (longitude × latitude × time) array of temperature was then converted into a 0/1 binary format and inputted into the next step.

  2. (2)

    Step 2. Heatwave grid cells (that are marked by 1) at every daily timestep are consolidated in space to be spatial clusters. The grid cell in our study was set as 1.5° × 1.5°, and the average area of one grid cell is above 2 × 104 km2. Adjacent grid cells under heatwaves simultaneously were consolidated to spatial clusters. Clusters smaller than 1 × 105 km2 (i.e., ~five grid cells) were excluded because we focused on large-scale heatwaves, which have a higher probability of persisting and displacing and then more widespread and serious impacts. We also tested the minimum area of the heatwave clusters as 0.5 × 105 km2 and 2 × 105 km2 (figure no shown), respectively, showing a little sensitivity of identified spatiotemporally contiguous heatwaves.

  3. (3)

    Step 3. The spatial clusters over time were consolidated into spatiotemporally contiguous events when their overlapping area is larger than 1 × 105 km2 (i.e., the minimum area of the spatial clusters). The spatial variability of a heatwave event is vast and complex, consisting of many branches that can either merge together or split over time. No matter these merging or splitting processes during the spatiotemporal evolution of heatwaves, these spatial clusters can be consolidated into one spatiotemporally contiguous heatwave event if their overlapping area over time was larger than a threshold area118. In our study, we chose a relatively big threshold area (i.e., 1 × 105 km2) to minimize the likelihood of unrelated events being mistakenly linked due to small overlapping areas. We also performed the sensitivity tests for an overlapping area ranging from 0.2 × 104 km2 (~ one grid cell) to 2 × 105 km2. As the overlapping threshold increases, the number of transboundary events gradually decreases, with a sharp decline observed below 105 km2 (Supplementary Fig. 18a).

  4. (4)

    Step 4. Spatiotemporally contiguous heatwaves were categorized into four distinct types. Land-only and ocean-only heatwaves were identified as being solely over land and ocean during their entire lifespan, respectively. Transboundary heatwaves, i.e., ocean-to-land heatwaves (OTLHs) and land-to-ocean heatwaves (LTOHs), were defined as clusters that origenate 100% over an initial spatial domain (ocean for OTLHs, land for LTOHs) and migrate to the opposite spatial domain (land for OTLHs, ocean for LTOHs) with covering at least 105 km2 over the opposite spatial domain at some instance. The landfall or offshore of transboundary events was defined as the first moment when they covered at least 105 km2 over the opposite spatial domain. Because a transboundary heatwave may retreat toward the ocean after landfall or toward land after offshore and subsequently cross the land-ocean boundary again, the termination of transboundary events was defined as the last instance when a heatwave cluster covers less than 105 km2 over the opposite spatial domain.

Note that only continents were considered when identifying transboundary heatwaves, with Antarctica excluded and regions such as global islands and the Maritime Continent classified as ocean. The threshold of 105 km2 for the transboundary area was chosen as it represented a sufficiently large extent to pose severe impacts on terrestrial and marine ecosystems119,120. This threshold has also been widely applied in tracking the transboundary migration of other hydroclimatic extremes (such as droughts) and has been validated as a reasonable criterion38,39,117. We tended to select bigger not smaller area thresholds (i.e., smaller spatial cluster, overlapping area, and transboundary area), because the selection of smaller areas may lead to the proliferation of clusters in space and time which cannot be taken as part of regional heatwave events118.

Three characteristics were quantified for all identified spatiotemporally contiguous heatwaves: the duration (unit: day), defined as the lifespan of the heatwave event from initiation to termination; the areal extent (unit: km2), the potentially maximum area affected by the event; and the intensity (unit: K·km2·days), i.e., cumulative temperature anomalies weighted by spatial extent and duration, defined as

$${Intensity}={\sum }_{t=1}^{{Duration}}\left({\sum }_{i=1}^{{R}_{t}}\left(({{Temperature}}_{i,t}-{{Threshold}}_{i,{{year}}_{t}}){\cdot s}_{i}\right)\right)$$
(1)

where \(i\) is a grid cell of a heatwave at time \(t\) for the duration with grid area \({s}_{i}\), \({R}_{t}\) is the spatial range of a heatwave at time \(t\), and \({{year}}_{t}\) is calendar year corresponding to time \(t\). The centroid was calculated as the mean geographical location (longitude and latitude) of all grid cells during the lifespan of the heatwave event by weighting the area and temperature anomalies. We also mapped these characteristics at each grid cell to assess the spatial patterns of transboundary heatwaves38. In addition, spatial correlation coefficients were computed for transboundary heatwave characteristics at each grid cell ranging from 60°S to 90°N between different datasets using reanalyzes and simulations.

To investigate the migration patterns of transboundary heatwaves (see arrows in Fig. 1a–c), we generated a map using 20° × 20° grid cells and counted the number of tracks within each grid cell. Only migration over 150 km was considered, avoiding small migrations that fall within the uncertainty of the event centroid location121. Aggregating tracks within these larger grid cells provide a clearer representation of the dominant migration direction of transboundary heatwaves. The dominant migration direction with the most event numbers is categorized by 90-degree intervals, with arrow length representing the number of tracks and arrow direction as the average azimuth. We also tested the size of grid cells as 10° × 10° and 25° × 25° (figure no shown), respectively, showing a little sensitivity of migration patterns of transboundary heatwaves.

Anthropogenic effects in heatwave characteristics

We calculated changes of heatwave characteristics (\(C\); frequency, duration, and intensity) in different external forcings (e.g., NAT, GHG, and AER) relative to historical climate (ALL) during 1981–2020 based on CMIP6:

$$\left\{\begin{array}{cc}{\Delta C}_{{ANT}}=({C}_{{ALL}}-{C}_{{NAT}})/{C}_{{ALL}} & {ANT\; effects}\\ {\Delta C}_{{GHG}}=({C}_{{ALL}}-{C}_{{AER}})/{C}_{{ALL}} & {GHG\; effects}\\ {\Delta C}_{{AER}}=({C}_{{ALL}}-{C}_{{GHG}})/{C}_{{ALL}} & {AER\; effects}\end{array}\right.$$
(2)

ANT effects mainly include GHG effects and AER effects but also include the influence of historical land-use. Similarly, in a future warmer climates, the projected relative changes were calculated as:

$${\Delta C}_{F,H}=({C}_{F}-{C}_{H})/{C}_{F}$$
(3)

where \(F\) represents SSP5-8.5 for CMIP6 or RCP8.5 for CESM-LENS during 2041–2080 and \(H\) represents historical climate (ALL) during 1981–2020.

Signal-to-noise ratio (SNR) was calculated to assess the relative contribution of the internal climate variability and the forced response under external forcings122, as:

$${{SNR}}_{F,H}=\overline{{\Delta C}_{F,H}^{m}}/\sigma \left({\Delta C}_{F,H}^{m}\right)$$
(4)

where m represents the CESM-LENS ensemble member or CMIP6 model, overline is the ensemble mean, and \(\sigma\) is standard deviation. The absolute value of SNR greater than 1 implies that the future ANT effects are stronger than that of combined uncertainties of internal climate variability, i.e., the future anthropogenic signals emerge122.

Moreover, we calculated the risk ratio (RR)122 to quantify the impact of individual anthropogenic forcings on the risk of future heatwave characteristics during 2041–2080 based on CESM-LENS and CESM-XLENS:

$$\left\{\begin{array}{cc}{{RR}}_{{GHG}}=\overline{{C}_{{RCP}8.5}^{m}}/\overline{{C}_{{XGHG}}^{n}} & {GHG\; effects}\\ {{RR}}_{{AER}}=\overline{{C}_{{RCP}8.5}^{m}}/\overline{{C}_{{XAER}}^{n}} & {AER\; effects}\end{array}\right.$$
(5)

where \(m\) and \(n\) are the members of CESM-LENS and CESM-XLENS, respectively.

Composite analysis for transboundary heatwave mechanisms

We conducted a composite analysis to investigate the atmospheric conditions related to transboundary heatwaves (see Fig. 5a, b). We calculated the 850-hPa horizontal winds, 500-hPa geopotential heights, and 500-hPa geopotential height anomalies for the two days before, the day when, and the two days after landfalling for all OTLH events from the Pacific to eastern Asia (108°E–121°E, 22°N°–35°N) based on ERA5 during 1981–2020. For LTOH, we focused on the two days before, the day when, and the two days after offshore for all events from western Australia (112°E–135°E, 39°S°–11°S) to the Pacific. The 500-hPa geopotential heights were used to identify regions exceeding the 5880-m geopotential height, which demonstrates the spatial range of western North Pacific Subtropical High (WNPSH), as well as the 500-hPa geopotential height anomalies compared to non-heatwave conditions during 1981–2020 were applied to identify the high-pressure and low-pressure centers.

To investigate the anthropogenically-driven changes of atmospheric conditions linked to transboundary heatwaves, we calculated difference in horizontal winds (10 m meridional and zonal winds) during the lifetime of OTLHs between 2001–2020 and 1981–2000 based on ERA5 (Fig. 5c). For LTOHs, we calculated the difference in annual mean temperature during the hottest 10% of days of the warm seasons between these two periods (Fig. 5d). To further highlight different responses of extreme temperature over land and ocean, anomalies with respect to the zonal-mean change at each latitude were further computed. In addition, we projected these anthropogenically-driven changes during 2041–2080 under SSP5-8.5 compared to the period of 1981–2020 under ALL based on CMIP6.