# Decrease in air-sea CO2 fluxes caused by persistent marine heatwaves

Alexandre Mignot, Karina von Schuckmann, Peter Landschützer, Florent Gasparin, Simon van Gennip, Coralie Perruche, Julien Lamouroux, Tristan Amm

► To cite this version:

Alexandre Mignot, Karina von Schuckmann, Peter Landschützer, Florent Gasparin, Simon van Gennip, et al.. Decrease in air-sea CO2 fluxes caused by persistent marine heatwaves. Nature Communications, 2022, 13 (1), pp.4300. 10.1038/s41467-022-31983-0. hal-04889847

HAL Id: hal-04889847

https://hal.science/hal-04889847v1

Submitted on 16 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Decrease in air-sea CO₂ fluxes caused by persistent marine heatwaves

Received: 8 January 2021

Accepted: 12 July 2022

Published online: 25 July 2022

<figure><img src="image_3.png" /></figure>

Check for updates

Alexandre Mignot<fnref n="1" />✉, Karina von Schuckmann<fnref n="1" />, Peter Landschützer², Florent Gasparin<fnref n="1" />, Simon van Gennip<fnref n="1" />, Coralie Perruche<fnref n="1" />, Julien Lamouroux<fnref n="1" /> & Tristan Amm<fnref n="1" />

Regional processes play a key role in the global carbon budget. Major ocean CO₂ uptake at mid-latitudes counteracts CO₂ release in the tropics, which is modulated by episodes of marine heatwaves. Yet, we lack essential knowledge on persistent marine heatwaves, and their effect on the CO₂ sensitive areas. Here we show, using a 1985–2017 joint analysis of reconstructions, ocean reanalysis and in situ and satellite data, that persistent marine heatwaves occur in major CO₂ uptake and release areas. Average air-sea CO₂ flux density changes from persistent marine heatwaves are strongest in the Pacific Ocean with a 40 ± 9% reduction in CO₂ release in the tropics linked to ENSO, and a reduction in CO₂ uptake of 29 ± 11% in the North Pacific over the study period. These results provide new insights into the interplay of extreme variability and a critical regulating ocean ecosystem service, and pave the way for future investigations on its evolution under climate change.

Extreme events associated with unusually high-water temperature are ubiquitous in the global ocean. They can last from weeks to years, span from local to interbasin scale, and can reach depths of several hundreds of meters<fnref n="1" />⁻³. These so-called marine heatwaves (MHWs) occur due to either coupled air-sea interactions⁴⁻⁷, ocean internal processes such as horizontal and/or vertical circulation changes⁸, and are sometimes linked to large climate modes such as the El Niño-Southern Oscillation (ENSO)⁵. Over the past 35 years, MHWs have become longer-lasting, more intense and more extensive⁹⁻¹¹ very likely due to long-term anthropogenic change<fnref n="1" />,⁹,¹⁰,¹²,¹³.

Intense and long-lasting MHWs have been reported at different locations in the global ocean³,⁵. These include the 2013/2015 Northeast Pacific 'warm blob'⁴,¹⁴, the 1997/1998 El Niño¹⁵, the 2015/2016 Tasman Sea¹⁶ or the 2012 Northwest Atlantic¹⁷,¹⁸. The duration of these major events ranges from several months up to 2 years, they are associated with dramatic increase in sea surface temperature, that can exceed 5 °C in anomalies¹⁹, and they can extend over large regions, reaching sometimes ~ 10 M km²⁵. Due to their extreme nature, MHWs, and in particular the intense, persistent ones, pose a fundamental challenge for societies as they have devastating impacts on marine ecosystem and their services<fnref n="1" />,¹¹,²⁰.

The ocean acts as a net sink for atmospheric CO₂ and is absorbing about a ¼ of CO₂ anthropogenic emissions²¹ (2.6 ± 0.6 PgC/yr over the 2009–2018 period), thereby mitigating global warming. CO₂ entering the ocean is then redistributed horizontally over large distances and into deep ocean layers where it is then stored for long time scales²²⁻²⁴. The magnitude and direction of the air-sea CO₂ flux density (FₐCO₂) vary widely in space and time, and depend on hydrographic conditions, the ocean circulation system, biological net production and air-sea interactions. As a result, major CO₂ uptake areas are located at mid-latitudes, whereas CO₂ release takes place predominantly in upwelling areas such as the tropical ocean²⁵⁻²⁷ (Fig. 1a).

Persistent MHWs linked to ENSO³ affect the Tropical Pacific CO₂ source region and lead to significant reduction in CO₂ outgassing²⁸⁻³². However, we lack essential knowledge about how major MHWs events affect other oceanic CO₂ sources and sink regions. Here we investigate the interplay and the impact of intense and long-lasting MHWs on the air-sea CO₂ flux density at the global scale. The study is built on a combined use of reconstructions from 1985 to 2017, direct measurements, remote sensing data and an ocean reanalysis. We first present the regions where particularly intense and long-lasting MHWs most frequently occur. We then quantify the impact of these extreme ocean events on oceanic CO₂ sink and sources areas. We further examine the

<footnote n="1">¹Mercator Océan International, Toulouse, France. ²Max Planck Institute for Meteorology, Hamburg, Germany. ³e-mail: mignot@mercator-ocean.fr</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1 | Interplay of PMHWs and oceanic carbon CO₂ source and sink areas.**</figcaption></figure>

a Mean 1985–2017 air-to-sea CO₂ flux density ($F_{CO2}$) derived from the Copernicus Marine Service (CMEMS) observation-based product (see methods section). Negative values indicate oceanic uptake (blue), while positive values indicate oceanic outgassing (red) of CO₂. The black continuous/dashed contours represent critical CO₂ sink/source regions, i.e the regions where the mean 1985–2017 $F_{CO2}$ is lower/greater than –1/1 molC/m²/yr as proposed by Takahashi et al.³⁵. The gray points represent satellite Sea Surface Temperature grid points that have experienced at least 3 PMHWs from 1985 to 2017 (see text for details). b Trimmed average percent $F_{CO2}$ anomalies during PMHWs derived from an ensemble of four observation-based products of $F_{CO2}$ (see section methods) in critical oceanic CO₂

sinks and sources (plain and dashed contours in Fig. 1a) that are impacted by PMHWs. The percent $F_{CO2}$ anomalies correspond to the monthly $F_{CO2}$ anomalies divided by the monthly $F_{CO2}$ climatological values (see section methods). Negative values correspond to a reduction in both a source or a sink region. The 95% confidence interval for each trimmed average percent $F_{CO2}$ anomalies are indicated in the Supplementary Table 1. We verified that all trimmed average percent $F_{CO2}$ anomalies were significantly different from 0 using a Yuen's trimmed mean test⁵⁹. The ensemble mean and standard deviation are given in black. An additional 12% uncertainty resulting from uncertain gas exchange⁶⁰ has been added to the ensemble spread. The calculation of the percent $F_{CO2}$ anomalies is detailed in the method section.

interaction between these extreme ocean events and one critical oceanic CO₂ sink region in the North Pacific Ocean. Finally, we discuss these results with existing knowledge on the mechanisms in the Tropical Pacific region to obtain a large-scale view of the prevailing mechanisms driving coupled changes between one of the regulating ocean ecosystem service and extreme variability.

## Results

### Persistent marine heatwaves occurrence and oceanic CO₂ source and sink areas

In the Tropical Pacific, intense and long-lasting MHWs (hereinafter denoted persistent marine heatwaves, PMHWs) have a strong impact on air-sea CO₂ fluxes³¹. We propose a specific new criteria to identify such PMHWs at the global scale based on the duration and the mean Sea Surface Temperature (SST) anomaly during a MHW. We first detect all MHWs that occurred from 1985 to 2017 by applying a standard MHW detection algorithm² to NOAA gridded SST data derived from AVHRR sensor³³,³⁴ (see method section). The detection algorithm provides several metrics that describe MHWs, including the duration and the mean SST anomaly. Using these two metrics, we define PMHWs

as MHWs whose duration and mean SST anomaly are greater than the 95th percentile of their global historical distribution, i.e., duration > 38 days and mean SST anomaly > 2.3 degrees Celsius. Finally, we focus on the regions where PMHWs have appeared several times over the past three decades, and that represent a recurring event effecting the ocean CO₂ sink, similar to El Niño events in the Tropical Pacific. To do so, we only consider the points where PMHWs have re-occurred at least three times during the 1985–2017 period (gray points in Fig. 1a)–which correspond to 25% of all grid points that have experienced at least one PMHW.

PMHWs most frequently occur in the largest oceanic CO₂ source and sink areas in the near-global ocean. The analysis is restricted to within 60°S and 60°N to exclude polar regions where periods of ice longer than 5 days obstruct the MHW detection algorithm (see method section). Critical sink regions are located at mid-latitudes in the Northern and Southern Hemispheres (solid contours in Fig. 1a), while ocean CO₂ outgassing predominantly occurs in upwelling regions such as the Tropical Pacific (dashed contours in Fig. 1a) and correspond to regions where climatological $F_{CO2}$ is lower/greater (sinks/sources) than –1/1 molC/m²/year (solid and dashed contours in Fig. 1a), as proposed
<table><caption>**Table 1 | Ranges in longitude and latitude, total area, area and normalized area impacted by PMHWs of the four regions considered in this study**</caption><thead><tr><th>critical CO<sub>2</sub> sink/source regions</th><th>Longitude range</th><th>Latitude range</th><th>Total Area (10<sup>6</sup> km<sup>2</sup>)</th><th>Area impacted by PMHWs (10<sup>6</sup> km<sup>2</sup>)</th><th>Area impacted by PMHWs/ total Area (%)</th></tr></thead><tbody><tr><td>North Pacific CO<sub>2</sub> sink</td><td>123.5°E-121.5°W</td><td>23.5°N-59.5°N</td><td>17.7</td><td>3.8</td><td>21</td></tr><tr><td>Tropical Pacific CO<sub>2</sub> source</td><td>173.5°E-73.5°W</td><td>17.5°S-5.5°N</td><td>18.4</td><td>2.9</td><td>16</td></tr><tr><td>North Atlantic CO<sub>2</sub> sink</td><td>74.5°W-7.5°E</td><td>30.5°N-65.5°N</td><td>11.0</td><td>1.1</td><td>10</td></tr><tr><td>Mid-high latitude Southern Oceans CO<sub>2</sub> sink</td><td>179.5°W-122.5 W and 85.5°W-180°E</td><td>20.5°S-56.5°S</td><td>56</td><td>2.8</td><td>5</td></tr></tbody></table>

<figure><img src="image_5.png" /></figure>

**Fig. 2 | Processes that lead to a reduction in the oceanic uptake of CO₂ in the North Pacific during PMHWs.** a 2009–2017 average F<sub>CO2</sub> anomalies (black dot) and its Taylor decomposition (vertical bars). The contribution of temperature, Dissolved Inorganic Carbon (DIC), Alkalinity (ALK), salinity, wind and atmospheric partial pressure of CO₂ to F<sub>CO2</sub> anomalies observed during PMHWs in the North Pacific CO₂ sink region for the 2009–2017 period were calculated using a first order Taylor expansion derived from the biogeochemical reanalysis (see methods section). The “total” bar corresponds to the sum of all contributing terms and corresponds to the Taylor approximation of F<sub>CO2</sub> anomalies (black dot). The good agreement between the two implies that F<sub>CO2</sub> anomalies are well approximated by

the Taylor decomposition. The error bars correspond to the 95% confidence interval. b Contribution of horizontal, vertical diffusion, vertical advection, dilution and concentration due to freshwater fluxes, air-sea CO₂ flux density, a residual term and biological activity to the rate of change (tendency or trend) of DIC anomalies during PMHWs in the North Pacific CO₂ sink region for the 2009–2017 period (see methods section). The vertical bars represent the slope from linearly regressing each forcing term to the DIC anomalies trend²⁸. A linear regression slope close to 1 indicates that a particular term produces in-phase anomalies of comparable magnitude. A slope near zero indicates that the term is not important in generating anomalies. The error bars correspond to 95% confidence intervals.

by Takahashi et al.³⁵. The range in longitude and latitude of the four ocean CO₂ sink/source regions are indicated in Table 1. Note that the mid-high latitude Southern Oceans CO₂ sink region includes the uptake area in the South Atlantic, Pacific, and Indian Oceans between 20.5° and 56.5°S, respectively. The climatological F<sub>CO2</sub> values are illustrated for the 1985–2017 period using the observation-based Copernicus Marine Environment Monitoring Service (hereafter denoted CMEMS) product³⁶ (see method section). Unexpectedly, we find that

the regions with a strong occurrence of PMHWs (gray points in Fig. 1a) are mainly located in the strongest oceanic sources and sinks areas in the near-global ocean, particularly in the Pacific Ocean (Fig. 1a and Table 1). In the North Pacific and in the Tropical Pacific, the area impacted by PMHWs covers 3.8 and 2.9 * 10⁶ km², respectively, which corresponds to 21% and 16% of the total area of these basins. Whereas, in the North Atlantic and in the mid-high latitude Southern Oceans, the
MECHANIMS LEADING TO A REDUCTION IN THE AIR-SEA CO₂ FLUX DENSITY IN THE PACIFIC DURING PERSISTENT MARINE HEATWAVES

<figure><img src="image_6.png" /><figcaption>**Fig. 3** | Schematic presentation of the mechanisms driving the reduction in the air-sea CO₂ flux density in (a) the North Pacific CO₂ sink and (b) the Tropical Pacific CO₂ source regions. Red color indicates the thermal effect on the air-sea</figcaption></figure>

CO₂ flux density, the blue color is linked to impacts related to circulations changes associated with PMHWs such as anomalous horizontal and vertical advection. The gray color represents the normal conditions. See text for more details.

area impacted by PMHWs covers only 10% and 5% of the total area of these basins.

The PMHWs impact on FCO2 is quantified using an ensemble of four observation-based products of FCO2 from 1985 to 2017 (see method section). We use a common approach using the ensemble spread as first order uncertainty estimate³⁷. We would like to note though, that regionally, the uncertainty might be larger resulting from the lack of direct measurements. Previous studies, however, show that the uncertainty as a result of data sparsity is in the order of the ensemble spread adopted here³⁸. Additionally data products are limited by 1×1 degree (2.5×2 degree for the JENA product) resolution by design and we are unable to formally quantify the role of the resolution on the uncertainty of our study. Amongst the largest oceanic CO2 sinks and sources where PMHWs most frequently occur, the North Pacific and the Tropical Pacific are the most impacted. During PMHWs events, both areas suffer from a significant reduction in the air-sea CO2 flux density with a change in uptake (29 ± 11%) and outgassing (40 ± 9%) respectively. In contrast, in the North Atlantic and the mid-high latitude Southern Oceans CO2 sinks, the impact of PMHWs on FCO2 is close to 0 and negligible over the study period. We verified that, in these two regions, the small value is a general feature and is not just a cancellation of large positive and negative anomalies. The impact of PMHWs on FCO2 appears thus to be the most important in the Tropical and North Pacific, which is of considerable concern given their contributions to the global ocean carbon cycle²⁵⁻²⁷. The Tropical (dashed contour) and North Pacific (solid contour) represent a source and a sink of 0.52 and -0.60 PgC yr⁻¹, respectively which corresponds to ~36 and 42% of the annual near-global net flux while covering roughly 5% of the near-global ocean area.

Persistent marine heatwaves and the North Pacific CO₂ sink

We use a state-of-the-art ocean biogeochemical reanalysis³⁹ (see supplementary information), validated against FCO₂ reconstructions (Fig. S1) and in situ observations from Biogeochemical (BGC)-Argo floats⁴⁰,⁴¹ (Fig. S2), to understand the interaction between PMHWs and FCO₂ in the North Pacific from 2009 to 2017. Note that, the calculations

relative to the 2009–2017 period are performed on those North Pacific regions where PMHWs have re-occurred several times since 1985, i.e. the SST grid points that have experienced at least 3 PMHWS from 1985 to 2017. The exchange of CO₂ between the ocean and the atmosphere is driven by six variables⁴²: wind, upper-ocean temperature, salinity, dissolved inorganic carbon (DIC) and alkalinity (ALK) as well as the atmospheric partial pressure of CO₂. Temperature, salinity, DIC and ALK are influenced by hydrographic conditions, the ocean circulation system, and air-sea interactions. In addition, DIC and ALK are also influenced by the biological net production⁴². As observation-based products of FCO₂ do not include these variables, the reanalysis becomes essential to pursue the analysis. The BGC reanalysis combines ocean circulation and biogeochemistry models together with in situ and satellite observations to provide a high degree of bio-physical realism⁴³. The reanalysis skill is validated against the ensemble of observation-based products over the 2009–2017 period and in situ observations from an array of BGC-Argo floats during the 2013/2015 ‘warm blob’ PMHW (see supplementary Information). The reanalysis shows good agreement with the observation-based products in estimating FCO₂ anomalies due to PMHWs in the North Pacific. The reanalysis also agrees well with the float observations in reproducing anomalies in the four oceanic drivers known to control FCO₂ (temperature, salinity, DIC and ALK) during the ‘warm blob’.

During PMHW events, the reduction in FCO₂ is the result of higher-than-usual temperature and negative DIC anomalies. We calculate a first-order Taylor series expansion of FCO₂ anomalies to determine the contribution of the four oceanic drivers²⁸,⁴⁴,⁴⁵ (see method section). The Taylor decomposition (Fig. 2a) reveals that the reduced uptake of CO₂ during PMHWs in the North Pacific mainly result from the contribution of temperature (1.43 ± 0.02 molC/m²/yr), DIC anomalies (-0.81 ± 0.01 molC/m²/yr) and to a lesser extent ALK anomalies (-0.23 ± 0.01 molC/m²/yr). The contribution from salinity, wind and the atmospheric partial pressure of CO₂ anomalies are small, and can be considered negligible. Sea surface warming during PMHWs reduces the solubility of CO₂ in the ocean resulting in a reduced uptake of CO₂. In contrast, the decrease in DIC associated with a small increase in ALK (Fig. S3) enhances the
uptake of CO₂ and as such offsets the thermal effect to the extent that the final FCO₂ anomaly is ~4 times smaller than it would have been from the thermal effect alone.

Next, we investigate the mechanisms leading to negative DIC anomalies. We examine the processes that drive the rate of change (tendency or trend) of DIC anomalies during PMHWs over the 2009–2017 period. The budget (or forcing) terms in the DIC trend equation consist of: horizontal and vertical advection, vertical diffusion, the air-sea CO₂ flux density, biological activity, dilution and concentration due to freshwater fluxes and a residual term (see method section). To highlight the contribution from each process to the DIC anomalies trend, we follow the method of Doney et al.²⁸ and examine the slope through linear regression of each forcing term to the DIC anomalies trends (Fig. 2b) (we verify that the intercept is approximately 0 because the average of the forcing term anomalies is null). A slope close to 1 indicates that a particular forcing term produces in-phase anomalies of comparable magnitude to DIC anomalies trend. In contrast, a slope near zero indicates that the term is not important, and a negative slope that the term produces out of phase anomalies.

Horizontal advection is the main driver for DIC anomalies. The linear regression slope of horizontal advection on DIC anomalies trends is the largest (0.70 ± 0.02 (unitless)) whereas the slopes of the other forcing terms are much smaller (<0.16 (unitless) for vertical diffusion anomalies and lower than 0.06 (unitless) for all the other terms). Furthermore and consistently with Ayers and Lozier⁴⁶ and Gruber et al.²³, the reanalysis shows that, on average, there is a net horizontal divergence of DIC in the North Pacific CO₂ sink region (data not shown). The reanalysis suggests that the lateral removal of DIC is further accentuated during PMHWs, causing a decrease in DIC. In the North Pacific, extreme MHWs are associated with changes in horizontal advection due to wind speed modification⁴,¹⁹. This suggests that similar processes could potentially both drive thermal and DIC changes during PMHWs in the North Pacific. Given the importance of the horizontal transport of DIC in mitigating the impact of PMHWs on the uptake of CO₂, we propose that studies should address how PMHWs and ocean circulation are interconnected in this region.

## Discussion

We show that PMHWs (> 38 days and > 2.3 °C anomalies as defined in this study) most frequently occur in oceanic regions of major importance for the global carbon cycle: the Tropical Pacific CO₂ source area, and the CO₂ sink regions of the North Pacific, the North Atlantic and the mid-high latitude Southern Ocean. However, over the study period 1985–2017, PMHWs have impacted air-sea CO₂ exchange only in the North and tropical Pacific CO₂ sensitive areas. The processes of this interplay are provided in the schematic of Fig. 3.

In the North Pacific CO₂ sink, PMHW events cause a reduction in the CO₂ sink as a result of the net effect of two competing mechanisms: extreme higher-than-average temperature and anomalous DIC advection. The former causes a reduction in the solubility of CO₂ in ocean water, thereby reducing the ocean uptake of CO₂ whereas the latter increases the ocean CO₂ uptake—through decreased levels of DIC driven by horizontal advection—and as such attenuates the impact of the thermal effect. Overall, the thermal effect dominates the advection effect in our study, leading to a net reduction in the air-to-sea CO₂ flux density during a PMHW event of about 29 ± 11% (Fig. 3a).

In the Tropical Pacific where CO₂ release takes place, the CO₂ outgassing is significantly attenuated during PMHWs with a reduction in CO₂ release from the ocean to the atmosphere of about 40 ± 9%. In this region, PMHWs are associated with ENSO⁵ and previous studies have investigated the mechanisms explaining this change, which is mainly driven by a change in vertical ocean circulation²⁸⁻³². During PMHWs (Fig. 3b), eastward propagating Kelvin waves that depress the thermocline in the east together with a

concurrent weakening of easterly winds, and the extension of the western Pacific warm pool to the east, reduce the upwelling of DIC leading to a net decrease in the sea-to-air CO₂ flux density. We have only considered the average sea-to-air CO₂ flux density response to PMHWs in the North and Tropical Pacific. Future studies should investigate the spatial variability of the response in these two regions. The results in the North Pacific CO₂ sink complete previous studies of ENSO-related PMHWs in the Tropical Pacific, and together with the new results obtained in this study thus provide a comprehensive view on the interplay between PMHWs and CO₂ sensitive areas in the North and Tropical Pacific as illustrated in Fig. 3.

Over the 1985–2017 period, considering the duration and affected areas in the North Pacific and Tropical Pacific, we derive an integrated flux anomaly linked to PMHWs of 13 ± 6 TgC (Teragrams of Carbon = 10¹² grams of Carbon) and –118 ± 25 TgC in the North Pacific and Tropical Pacific respectively. The climatological integrated fluxes expected during the same time and within the same area are –119 ± 12 TgC and 307 ± 31 TgC respectively. This corresponds to a flux anomaly of 11 ± 5% and 39 ± 11% for the respective regions. While the effect in absolute terms appears small at first sight, particularly in the North Pacific (only about 0.5% of all annual marine uptake of anthropogenic CO₂³⁷), increasing duration and intensity of the heat waves become relevant for closing the current carbon budget imbalance (currently in the order of ±100–500 TgC/yr for the past decades and recent years (see Table 6 in Friedlingstein et al.³⁷) and should therefore receive attention in future budgets.

In this study, we have introduced the concept of intense and long-lasting MHWs, here defined as PMHWs, with a criterion based on the duration and the mean SST anomaly (i.e., duration > 38 days and mean SST anomaly > 2.3 degrees Celsius). Our results have shown that PMHWs occur predominantly in the major CO₂ source and sink areas. Shorter-lasting MHW events are more widespread than PMHWs¹⁰,¹⁹ and further ensemble-based analyses are needed to understand on how their impact might affect air-sea CO₂ exchanges at global scale. Except for one single product (JENA), all other observational-based products are available at monthly resolution which makes such investigation not possible. Future analyses based on an ensemble of CO₂ products of high temporal resolution (i.e. < 30 day) would be needed to study the integrated impact of shorter extreme events on ocean CO₂ fluxes.

In the North Atlantic and the mid-high latitude Southern Ocean CO₂ sinks, no change is observed for the CO₂ fluxes due to PMHWs, suggesting that the thermal effect and non-thermal effects cancel each other out. To further analyze this hypothesis, the use of reanalyses is needed. Currently, long-term in-situ observations for pCO₂, temperature, salinity, DIC and ALK are lacking in these areas preventing the establishment of robust air-sea CO₂ flux products, climatologies and reanalysis skill validation for such a process study. Hence, PMHW-induced impacts on CO₂ fluxes mechanisms in the Atlantic and mid-high latitude Southern Ocean remain unanswered, and increased monitoring efforts are needed accordingly.

There has been a statistically significant increase (P=.08) in the attenuation in FCO₂ due to PMHWs between the 1985–1995 and the 2007–2017 period in the North Pacific, and we have compared the distribution of the ensemble of the 4 observation-based products. On the contrary, there has been no statistically significant change in the Tropical Pacific (Fig. 4). Similarly, PMHWs have also increased in intensity in the North Pacific over the last decades (linear trend= 0.11 ± 0.03 °C/decade, P=.002), while their strength remains similar in the Tropics (linear trend= 0.04 ± 0.09 °C/decade, P=.66) (Fig. S4). Based on the results of our process study, we can develop the hypothesis that the reported increase in the intensity of PMHWs has potentially amplified the outgassing of CO₂ over the 1985–2017 period, and that the competing mechanism, i.e., anomalous advection of DIC, was unable to counter-interact the thermal effect over this time scale. However, we cannot test such hypothesis as it would require a
<figure><img src="image_7.png" /><figcaption>Fig. 4 | Evolution of FCO2 anomalies due to PMHWs during the 1985–2017 period. Trimmed average percent FCO2 anomalies during PMHWs for three time-periods (1985–1995, 1996–2006, 2007–2017) derived from an ensemble of 4 observation-based products of FCO2 (see section methods) in a North Pacific CO2</figcaption></figure>

sink and b the Tropical Pacific CO2 source regions. The ensemble mean and standard deviation are given in black. An additional 12% uncertainty resulting from uncertain gas exchange⁶⁰ has been added to the ensemble spread. The calculation of the percent FCO2 anomalies is detailed in the method section.

decomposition of FCO2 and DIC budgets from 1985 to 2017, the latter being currently not estimated by the reanalysis over this period. MHWs are projected to become stronger, more frequent and longer lasting in a warming climate¹,⁹,¹²,¹³. Therefore, it is crucial to understand how PMHWs and FCO2 interact over longer time scales if we want to further unravel the evolution of the oceanic carbon cycle under climate change.

## Methods

### Observation-based products of FCO2

In this study, we use an ensemble of four observation-based products to quantify the impact of PMHWs on FCO2 in the three largest oceanic CO2 sink and the largest oceanic CO2 source in the Tropical Pacific from 1985 to 2017. Here, we provide a brief outline of the chosen products. More detail can be found in their respective publications.

The first observation-based product, from the Max Plank Institute for Meteorology (hereinafter denoted MPI)⁴⁷,⁴⁸, is based on a self-organizing map-feed-forward network that reconstructs the sea surface partial pressure of CO2 (spCO2) from various environmental predictor data. In a first step, the ocean is divided into biogeochemical regions of similar spCO2 properties (making use of a spCO2 climatology) and in a second step the non-linear relationship between auxiliary driver data and sparse observations is reconstructed to fill measurement gaps. The period of analysis is from 1982 to 2019 at monthly intervals and with a spatial resolution of 1° × 1°. It is based on a collection of ship and mooring spCO2 measurements assembled by the Surface Ocean CO2 Atlas (SOCAT) version 2020⁴⁹⁻⁵².

The second observation-based product, from Copernicus Marine Environmental Monitoring Service (hereinafter denoted CMEMS)³⁶, is from an ensemble-based forward feed neural network that reconstruct change in spCO2 from environmental predictor data. The period of analysis is from 1985 to 2018 at monthly intervals and with a spatial resolution of 1° × 1°. It is based on a collection of ship and mooring spCO2 measurements assembled by the Surface Ocean CO2 Atlas (SOCAT) version 2019⁴⁹⁻⁵².

The third observation-based product, from the Council for Scientific and Industrial Research (hereinafter denoted CSIR)⁵³, is from a machine-learning ensemble average of six two-step clustering-regression models that reconstruct change in spCO2 from environmental predictor data. The period of analysis is from 1982 to 2019 at monthly intervals and with a spatial resolution of 1° × 1°. It is based on a collection of ship and mooring spCO2 measurements assembled by the Surface Ocean CO2 Atlas (SOCAT) version 2019⁴⁹⁻⁵².

The fourth observation-based product, from the Max Plank Institute for Biogeochemistry (hereinafter denoted JENA)⁵⁴, is from an observation-driven ocean mixed-layer scheme that reconstruct change in spCO2 by fitting a data-driven diagnostic model of ocean mixed-layer biogeochemistry to surface-ocean CO2 partial pressure data from the SOCAT version 2019⁴⁹⁻⁵². The period of analysis is from 1957 to 2019 at daily intervals and with a spatial resolution of 2.5° × 2°. The daily fields were averaged into monthly fields.

Finally, to evaluate the skill of the BGC reanalysis in estimating FCO2 anomalies associated with PMHWs in the North Pacific CO2 sink, we use an additional observation-based product, from the Japan Meteorology Agency (hereinafter denoted JMA)⁵⁵, which is excluded from the spCO2 ensemble as its period of analysis is shorter than the previously listed products, i.e. from 1990 to 2018. This product is based on multiple linear regressions that reconstruct change in spCO2 from a set of environmental drivers. The temporal resolution is monthly intervals and with a spatial resolution of 1° × 1°. It is based on a collection of ship and mooring spCO2 measurements assembled by the Surface Ocean CO2 Atlas (SOCAT) version 2019⁴⁹⁻⁵².

### Estimates of the air-to-sea flux density of CO2 from spCO2 data

In the five observation-based products and the BGC reanalysis, the air-to-sea CO2 flux density (FCO2) is generated from spCO2 data using the gas exchange formulation³⁶,

$$ F_{CO2} = k\alpha(spCO2 - pCO2_{atm}), \quad (1) $$

where α is the CO2 solubility in seawater, k, a gas transfer coefficient, pCO2atm is the atmospheric partial pressure of CO2 and spCO2 is the sea surface partial pressure of CO2. Here, negative values of FCO2 indicate uptake of CO2 from the atmosphere to the ocean, while positive values indicate outgassing of CO2 from the ocean to the atmosphere. Each product performs their own calculation of the fluxes and the methods are described in the respective publications.

### Calculation of 2009–2017 monthly anomalies

In the reanalysis and the observation-based products, monthly anomalies (hereinafter denoted with a prime) are computed by removing a climatological value (hereinafter denoted with an overbar). The climatological value corresponds to the sum of a long-term linear trend and a monthly mean value. The monthly mean values are computed from the detrended monthly data.
## Calculation of 1985–2017 percent FCO2 anomalies

The percent FCO2 anomalies during PMHWs and for the 1985–2017 period correspond to the monthly FCO2 anomalies divided by the monthly FCO2 climatological values. The anomalies and climatological values were computed following the method detailed previously, with the exception that monthly mean values were only computed from the detrended monthly data over the 1985–1995 period. During this decade, the number of PMHWs per year, and near globally, were the lowest of the 1985–2017 period (see Fig. S5). By calculating the anomalies relative to this “reference” decade, we make sure that the percent FCO2 anomalies represent a change with respect to oceanic conditions not impacted by PMHWs.

In Fig. 1b, we represent, in each CO2 sink/source, an ensemble of four 1985–2017 trimmed mean percent FCO2 anomalies derived from the observation-based products. We use the trimmed mean instead of the mean because it is a robust estimator of central tendency and provides a better estimation of the location of the bulk of the data than the mean when the distribution is asymmetric, which is the case here. More precisely, we use a 5% trimmed mean, i.e., the lowest 5% and the highest 5% of the data are excluded.

## Taylor expansion of FCO2 anomalies

To determine the driving mechanisms causing FCO2 anomalies during PMHWs in the North Pacific CO2 sink, we calculate a first-order Taylor series expansion of FCO2 anomalies in terms of its driving parameters (i.e., wind, upper ocean temperature, salinity, dissolved inorganic carbon (DIC), alkalinity (ALK) and the atmospheric partial pressure of CO2)28,44.

First, we performed the linear Taylor decomposition of Eq. (1):

$$F_{CO_2} \approx (k\alpha)'(spCO_2 - pCO_{2atm}) + (\overline{k\alpha})spCO_2' - (\overline{k\alpha})pCO_{2atm}'. \quad (2)$$

The right-hand-side terms represent the contribution to FCO2' of gas transfer and solubility anomalies, atmospheric pCO2 anomalies and spCO2 anomalies. Note that the temperature dependence of k and α cancel each other, and (kα)' is mainly driven by variations in wind speed28,44.

The spCO2 anomalies are further decomposed into contributions from sea surface temperature anomalies (SST'), sea surface dissolved inorganic carbon anomalies (SDIC'), sea surface alkalinity anomalies (SALK') and sea surface salinity anomalies (SSS'), neglecting the second-order terms28,44,45,57:

$$spCO_2' \approx \frac{\partial spCO_2}{\partial SDIC} SDIC' + \frac{\partial spCO_2}{\partial SALK} SALK' + \frac{\partial spCO_2}{\partial SST} SST' + \frac{\partial spCO_2}{\partial SSS} SSS'. \quad (3)$$

Substituting Eq. (3) into Eq. (2) gives the contributions of all parameters to FCO2' in a single expression:

$$F_{CO_2}' = (k\alpha)'(spCO_2 - pCO_{2atm}) + (\overline{k\alpha})\left(\frac{\partial spCO_2}{\partial SDIC}SDIC' + \frac{\partial spCO_2}{\partial SALK}SALK' + \frac{\partial spCO_2}{\partial SST}SST' + \frac{\partial spCO_2}{\partial SSS}SSS'\right) - (\overline{k\alpha})pCO_{2atm}'. \quad (4)$$

Here, the sea surface quantities correspond to the quantities at the first level of the ocean reanalysis estimates (z ~ -0.50 m). Following Doney et al.28, the partial derivatives in Eq. (4) were computed off-line at each grid point, taking SDIC as an example, as:

$$\frac{\partial spCO_2}{\partial SDIC} \approx \frac{spCO_2(SDIC, SALK, SST, SSS) - spCO_2(SDIC, SALK, SST, SSS)}{SDIC'} \quad (5)$$

where spCO2 values are calculated using the seacarb program for R (https://CRAN.R-project.org/package=seacarb).

## DIC anomalies budget

In our study, we show that DIC anomalies play a significant role in controlling FCO2 anomalies during PMHWs. We therefore conduct a DIC anomalies budget to elucidate what processes controlled DIC anomalies.

In the ocean reanalysis, the changes in DIC concentration with time are described by the following equation:

$$\frac{\partial DIC}{\partial t} = ADV_H + ADV_z + DIFF_z + SBC + F_{CO2} + B + r \quad (6)$$

where ADVH and ADVz are the horizontal and vertical advection of DIC respectively, DIFFz is the vertical diffusion of DIC, SBC are the freshwater fluxes that dilute or concentrate DIC, FCO2 is the air-sea CO2 flux density, B is the biological activity that consumes or releases DIC (see details in Aumont et al.58), and r is the climatological damping (see supplementary information). Positive values result in a net increase in DIC. All terms were computed online on a daily basis and stored for monthly averages. The DIC tendency (rate of change or trend) equation (Eq. 6) is expressed as a function of monthly anomalies and averaged over the average mixing layer observed during PMHWs in the reanalysis (indicated by angle brackets), i.e. from the surface to h ~ 47 m:

$$\frac{\partial <DIC'>}{\partial t} = <ADV_H> + <ADV_z> + <DIFF_z> + <SBC> + \frac{<F_{CO2}>}{h} + <B> + <r>. \quad (7)$$

## Satellite sea surface temperature and marine heatwaves detection

MHWs locations, dates of onset and durations were derived from the global daily remotely sensed National Ocean Atmospheric Administration (NOAA) Optimum Interpolation sea surface temperature V2, 1/4° gridded data over 1982–201733,34. This dataset is derived from the advanced very high-resolution radiometer (AVHRR).

We apply a standard MHW detection algorithm2 to the gridded SST data. More specifically, a warm event is considered as a MHW if it lasts for 5 or more days, with sea surface temperatures warmer than the 90th percentile based on a 1983–2012 historical climatology. The MHW detection algorithm is usually not performed on grid cells with periods of ice coverage longer than 5 days10. We therefore restrict our analysis to the area between 60°S and 60°N. For each MHW detected, the date of onset, duration and mean sea surface temperature anomaly are estimated by the MHW detection algorithm.

## Calculation of anomalies associated with PMHWs

For each PMWH detected, the monthly anomalies were extracted at the model or observation-based products grid-point the closest to the PMHW location and for the entire duration of the PMHW. Then, to match the temporal resolution of the PMHW, the extracted anomalies were resampled from monthly to daily frequency through linear interpolation. The interpolated values were then averaged over the duration of the PMHW to give a single value, consistently with the other metrics derived from the MHW detection algorithm.

## Data availability

The reanalysis data can be downloaded from the Copernicus Marine Environmental Monitoring Service (https://resources.marine.copernicus.eu/product-detail/GLOBALengerAnalysisForecastBIO_001_028/information). The DIC budget terms data are available upon request from the corresponding author. The BGC-Argo data were
downloaded from the Argo Global Data Assembly Centre in France (ftp://ftp.ifremer.fr/argo/). The observation-based product are available from the Surface Ocean pCO2 Mapping Intercomparison website (http://www.bgc-jena.mpg.de/SOCOM/). The SST data are provided by NOAA/ESR/PSL at https://psl.noaa.gov/data/gridded/data.noaa.oisst.v2.highres.html.

## References

1. Collins, M. et al. Extremes, abrupt changes and managing risk. In *IPCC Special Report on the Ocean and Cryosphere in a Changing Climate* (2019).

2. Hobday, A. J. et al. A hierarchical approach to defining marine heatwaves. *Prog. Oceanogr.* **141**, 227–238 (2016).

3. Holbrook, N. J. et al. Keeping pace with marine heatwaves. *Nat. Rev. Earth Environ.* **1**, 482–493 (2020).

4. Bond, N. A., Cronin, M. F., Freeland, H. & Mantua, N. Causes and impacts of the 2014 warm anomaly in the NE Pacific. *Geophys. Res. Lett.* **42**, 3414–3420 (2015).

5. Holbrook, N. J. et al. A global assessment of marine heatwaves and their drivers. *Nat. Commun.* **10**, 1–13 (2019).

6. Lee, T. et al. Record warming in the South Pacific and western Antarctica associated with the strong central-Pacific El Niño in 2009–10. *Geophys. Res. Lett.* **37**, L19704 (2010).

7. Rodrigues, R. R., Taschetto, A. S., Sen Gupta, A. & Foltz, G. R. Common cause for severe droughts in South America and marine heatwaves in the South Atlantic. *Nat. Geosci.* **12**, 620–626 (2019).

8. Behrens, E., Fernandez, D. & Sutton, P. Meridional oceanic heat transport influences marine heatwaves in the Tasman Sea on interannual to decadal timescales. *Front. Mar. Sci.* **6**, 228 (2019).

9. Frölicher, T. L., Fischer, E. M. & Gruber, N. Marine heatwaves under global warming. *Nature* **560**, 360–364 (2018).

10. Oliver, E. C. J. et al. Longer and more frequent marine heatwaves over the past century. *Nat. Commun.* **9**, 1324 (2018).

11. Smale, D. A. et al. Marine heatwaves threaten global biodiversity and the provision of ecosystem services. *Nat. Clim. Chang.* **9**, 306–312 (2019).

12. Laufkötter, C., Zscheischler, J. & Frölicher, T. L. High-impact marine heatwaves attributable to human-induced global warming. *Science* **369**, 1621–1625 (2020).

13. Oliver, E. C. J. Mean warming not variability drives marine heatwave trends. *Clim. Dyn.* **53**, 1653–1659 (2019).

14. Di Lorenzo, E. & Mantua, N. Multi-year persistence of the 2014/15 North Pacific marine heatwave. *Nat. Clim. Change* **6**, 1042–1047 (2016).

15. McPhaden, M. J. Genesis and evolution of the 1997–98 El Nino. *Science* **283**, 950–954 (1999).

16. Oliver, E. C. J. et al. The unprecedented 2015/16 Tasman Sea marine heatwave. *Nat. Commun.* **8**, 16101 (2017).

17. Chen, K., Gawarkiewicz, G. G., Lentz, S. J. & Bane, J. M. Diagnosing the warming of the Northeastern U.S. Coastal Ocean in 2012: A linkage between the atmospheric jet stream variability and ocean response. *J. Geophys. Res. Oceans* **119**, 218–227 (2014).

18. Mills, K. et al. Fisheries management in a changing climate: lessons from the 2012 ocean heat wave in the Northwest Atlantic. *Oceanog* **26**, 191–195 (2013).

19. Sen Gupta, A. et al. Drivers and impacts of the most extreme marine heatwave events. *Sci. Rep.* **10**, 19359 (2020).

20. Ummenhofer, C. C. & Meehl, G. A. Extreme weather and climate events with ecological relevance: a review. *Philos. Trans. R. Soc. B* **372**, 20160135 (2017).

21. Friedlingstein, P. et al. Global carbon budget 2019. *Earth Syst. Sci. Data* **11**, 1783–1838 (2019).

22. Sabine, C. L. The oceanic sink for anthropogenic CO2. *Science* **305**, 367–371 (2004).

23. Gruber, N. et al. The oceanic sink for anthropogenic CO2 from 1994 to 2007. *Science* **363**, 1193–1199 (2019).

24. Ciais, P. et al. Carbon and other biogeochemical cycles. In *Climate change 2013: the physical science basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change* (Cambridge University Press, 2013).

25. Gruber, N. et al. Oceanic sources, sinks, and transport of atmospheric CO2. *Global Biogeochem. Cycles* **23**, GB1005 (2009).

26. Landschützer, P., Gruber, N., Bakker, D. C. E. & Schuster, U. Recent variability of the global ocean carbon sink. *Glob. Biogeochem. Cycles* **28**, 927–949 (2014).

27. Takahashi, T. et al. Climatological mean and decadal change in surface ocean pCO2, and net sea-air CO2 flux over the global oceans. *Deep Sea Res. Part II: Topical Stud. Oceanogr.* **56**, 554–577 (2009).

28. Doney, S. C. et al. Mechanisms governing interannual variability in upper-ocean inorganic carbon system and air–sea CO2 fluxes: Physical climate and atmospheric dust. *Deep Sea Res. Part II: Topical Stud. Oceanogr.* **56**, 640–655 (2009).

29. Feely, R. A. et al. Seasonal and interannual variability of CO2 in the equatorial Pacific. *Deep Sea Res. Part II: Topical Stud. Oceanogr.* **49**, 2443–2469 (2002).

30. Le Quéré, C., Orr, J. C., Monfray, P., Aumont, O. & Madec, G. Interannual variability of the oceanic sink of CO2 from 1979 through 1997. *Glob. Biogeochemical Cycles* **14**, 1247–1265 (2000).

31. Liao, E., Resplandy, L., Liu, J. & Bowman, K. Amplification of the ocean carbon sink during El Niños: role of poleward Ekman transport and influence on atmospheric CO₂. *Global Biogeochem. Cycles* https://doi.org/10.1029/2020GB006574 (2020).

32. McKinley, G. A., Follows, M. J. & Marshall, J. Mechanisms of air-sea CO2 flux variability in the equatorial Pacific and the North Atlantic. *Global Biogeochem. Cycles* **18**, GB2011 (2004).

33. Banzon, V., Smith, T. M., Chin, T. M., Liu, C. & Hankins, W. A. long-term record of blended satellite and in situ sea-surface temperature for climate monitoring, modeling and environmental studies. *Earth Syst. Sci. Data*, **8**, 165–176 (2016).

34. Reynolds, R. W. et al. Daily high-resolution-blended analyses for sea surface temperature. *J. Clim.* **20**, 5473–5496 (2007).

35. Takahashi, T. et al. Global sea-air CO2 flux based on climatological surface ocean pCO(2), and seasonal biological and temperature effects. *Deep-Sea Res. Part II-Top. Stud. Oceanogr.* **49**, 1601–1622 (2002).

36. Denvil-Sommer, A., Gehlen, M., Vrac, M. & Mejia, C. LSCE-FFNN-v1: a two-step neural network model for the reconstruction of surface ocean pCO2 over the global ocean. *Geosci. Model Dev.* **12**, 2091–2105 (2019).

37. Friedlingstein, P. et al. Global Carbon Budget 2020. *Earth Syst. Sci. Data* **12**, 3269–3340 (2020).

38. Gloege, L. et al. Quantifying Errors in Observationally Based Estimates of Ocean Carbon Sink Variability. *Global Biogeochem. Cycles* **35**, e2020GB006788 (2021).

39. CMEMS. https://resources.marine.copernicus.eu/product-detail/GLOBALenger_sciences_001_028/information

40. Claustre, H., Johnson, K. S. & Takeshita, Y. Observing the global ocean with biogeochemical-argo. *Annu. Rev. Marine Sci.* **12**, 23–48 (2020).

41. Roemmich, D. et al. On the future of Argo: a global, full-depth, multi-disciplinary array. *Front. Mar. Sci.* **6**, 439 (2019).

42. Williams, R. G. & Follows, M. J. Ocean dynamics and the carbon cycle: Principles and mechanisms. (Cambridge University Press, 2011).

43. Fennel, K. et al. Advancing marine biogeochemical and ecosystem reanalyses and forecasts as tools for monitoring and managing ecosystem health. *Front. Mar. Sci.* **6**, 89 (2019).
44. Lovenduski, N. S., Gruber, N., Doney, S. C. & Lima, I. D. Enhanced CO2 outgassing in the Southern Ocean from a positive phase of the Southern Annular Mode. *Glob. Biogeochem. Cycles* **21**, n/a–n/a (2007).

45. Brady, R. X., Lovenduski, N. S., Alexander, M. A., Jacox, M. & Gruber, N. On the role of climate modes in modulating the air–sea CO2 fluxes in eastern boundary upwelling systems. *Biogeosciences* **16**, 329–346 (2019).

46. Ayers, J. M. & Lozier, M. S. Unraveling dynamical controls on the North Pacific carbon sink. *J. Geophys. Res. Oceans* **117**, CO1017 (2012).

47. Landschützer, P., Gruber, N. & Bakker, D. C. Decadal variations and trends of the global ocean carbon sink. *Glob. Biogeochemical Cycles* **30**, 1396–1417 (2016).

48. Landschützer, P., Gruber, N. & Bakker, D. C. E. An observation-based global monthly gridded sea surface pCO2 product from 1982 onward and its monthly climatology (NCEI Accession 0160558). https://doi.org/10.7289/V5Z899N6 (2017).

49. Bakker, D. C. E. et al. A multi-decade record of high-quality CO2 data in version 3 of the Surface Ocean CO2 Atlas (SOCAT). *Earth Syst. Sci. Data* **8**, 383–413 (2016).

50. Sabine, C. L. et al. Surface Ocean CO2 Atlas (SOCAT) gridded data products. *Earth Syst. Sci. Data Discuss.* **5**, 781–804 (2012).

51. Bakker, D. C. E. et al. An update to the Surface Ocean CO2 Atlas (SOCAT version 2). *Earth Syst. Sci. Data* **6**, 69–90 (2014).

52. Pfeil, B. et al. A uniform, quality controlled Surface Ocean CO2 Atlas (SOCAT). *Earth Syst. Sci. Data* **5**, 125–143 (2013).

53. Gregor, L., Lebehot, A. D., Kok, S. & Scheel Monteiro, P. M. A comparative assessment of the uncertainties of global surface ocean CO2 estimates using a machine-learning ensemble (CSIR-ML6 version 2019a) – have we hit the wall? *Geosci. Model Dev.* **12**, 5113–5136 (2019).

54. Rödenbeck, C. et al. Global surface-ocean pCO2 and sea-air CO2 flux variability from an observation-driven ocean mixed-layer scheme. *Ocean Sci.* **9**, 193–216 (2013).

55. Iida, Y. et al. Trends in pCO2 and sea-air CO2 flux over the global open oceans for the last two decades. *J. Oceanogr.* **71**, 637–661 (2015).

56. Wanninkhof, R. Relationship between wind speed and gas exchange over the ocean revisited: Gas exchange and wind speed over the ocean. *Limnol. Oceanogr. Methods* **12**, 351–362 (2014).

57. Turi, G., Lachkar, Z. & Gruber, N. Spatiotemporal variability and drivers of pCO2 and air–sea CO2 fluxes in the California Current System: an eddy-resolving modeling study. *Biogeosciences* **11**, 671–690 (2014).

58. Aumont, O., Ethé, C., Tagliabue, A., Bopp, L. & Gehlen, M. PISCES-v2: an ocean biogeochemical model for carbon and ecosystem studies. *Geosci. Model Dev.* **8**, 2465–2513 (2015).

59. Yuen, K. K. The two-sample trimmed t for unequal population variances. *Biometrika* **61**, 165–170 (1974).

60. Roobaert, A., Laruelle, G. G., Landschützer, P. & Regnier, P. Uncertainty in the global oceanic CO₂ uptake induced by wind

forcing: quantification and spatial analysis. *Biogeosciences* **15**, 1701–1720 (2018).

## Acknowledgements

This study has been conducted using the Copernicus Marine Service products. The BGC-Argo data were collected and made freely available by the International Argo program and the national programs that contribute to it (https://www.argo.jcommops.org). The Argo program is part of the Global Ocean Observing System.

## Author contributions

A.M., F.G., K.v.S. designed the study. A.M. and K.v.S. wrote the draft of the manuscript, together with F.G., P.L., C.P. and S.v.G. T.A. and S.v.G. ran the algorithm for the detection of MHWS. J.L. and C.P. ran the analysis. A.M. analyzed the data. All authors read and approved the final draft.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-022-31983-0.

**Correspondence** and requests for materials should be addressed to Alexandre Mignot.

**Peer review information** *Nature Communications* thanks Maria Calleja, Alex Sen Gupta and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher’s note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022