The ocean losing its breath under the heatwaves

Domain: earth_science
Published: 2024-08-09
Editorial outcome: minor_revision

Abstract
The world’s oceans are under threat from the prevalence of heatwaves caused by climate change. Despite this, there is a lack of understanding regarding their impact on seawater oxygen levels - a crucial element in sustaining biological survival. Here, we find that heatwaves can trigger low-oxygen extreme events, thereby amplifying the signal of deoxygenation. By utilizing in situ observations and state-of-the-art climate model simulations, we provide a global assessment of the relationship between the two types of extreme events in the surface ocean (0–10 m). Our results show compelling evidence of a remarkable surge in the co-occurrence of marine heatwaves and low-oxygen extreme events. Hotspots of these concurrent stressors are identified in the study, indicating that this intensification is more pronounced in high-biomass regions than in those with relatively low biomass. The rise in the compound events is primarily attributable to long-term warming primarily induced by anthropogenic forcing, in tandem with natural internal variability modulating their spatial distribution. Our findings suggest the ocean is losing its breath under the influence of heatwaves, potentially experiencing more severe damage than previously anticipated.

Introduction
The global ocean has been warming and losing its oxygen primarily as a consequence of human-induced changes to the Earth’s climate1–4. Alongside long-term trend of increasing temperature and declining oxygen concentration, risks of extreme condition in these properties, i.e., marine heatwave and low-oxygen extreme events, has emerged in the past few decades5,6. The substantial increases of these extreme signals in climate change could severely affect marine biota and ecosystems, which is recognized as a serious threat to sustainable future for the ocean and related human communities7–10.

Currently, an increasing number of studies investigated the impact of global warming on declining oxygen concentrations (i.e., ocean deoxygenation) over decadal and multidecadal periods11,12. The long-term warming is anticipated to reduce the solubility of oxygen and enhance biological consumption, and the warming-induced strengthening of stratification would lead to reduced ocean ventilation and a subsequent decrease in oxygen supply to the ocean interior13,14. There is also a rapidly burgeoning literature focusing on extreme warming events in the ocean, defined as marine heatwaves, which could cause widespread biological and socio-economic impacts15–19. Research on low-oxygen extreme events started relatively late, but these events have received increasing attention in recent years. Characterized by the abrupt oxygen decline in the ocean, the low-oxygen extreme events could cause species to reach their critical oxygen limits decades earlier than expected from the deoxygenation20, eventually resulting in the loss of biodiversity and habitats8. The emergence of these low-oxygen extreme events is influenced by multiple factors, including anomalous water mass distributions and ocean eddies etc.5,21. Importantly, case studies have suggested that marine heatwaves may act as catalysts for extreme low-oxygen events: the abnormally low dissolved oxygen levels were found during a strong marine heatwave known as “the Blob”22; In SW Atlantic, a marine heatwave with sea surface temperatures over 30 °C coincided with the first regional record of hypoxia23. Observations on the Baltic Sea reveal heatwaves can contribute to seawater oxygen deficiency, especially in the shallow areas24. These researches highlight the importance of exploring the coupling between heatwaves and extreme low-oxygen events. Many marine ecosystems are highly vulnerable to multiple stressors and as a consequence may suffer more severe damage if low-oxygen extreme and heatwave occur simultaneously or in close spatiotemporal proximity, as such joint occurrence has the potential to act synergistically, leading to mutually reinforcing impacts on the marine organisms25,26. The combination of high water temperature and low oxygen concentration has exerted influences on marine life such as fishes and scallops27–29. Nevertheless, so far we still lack a detailed knowledge about how heatwave and low-oxygen extreme coupling across the ocean. A global assessment of the modification of the compound events, as well as the drivers behind them, is urgently needed under ongoing climate change.

To fill this gap, here we conduct a comprehensive analysis of the relationship between marine heatwaves and low-oxygen extreme events in the global ocean, using in situ observations and a fully coupled Earth system model30 (Data descriptions and validations of model could be found in “Methods”). The evolutions of the compound events of the two extreme events are investigated in our study and a framework is established to identify the drivers contributing to their occurrences. We focus on extreme events within the uppermost ocean layer (0–10 m) because conditions in the surface ocean, such as temperature and oxygen levels, exert profound influences on marine ecosystems, encompassing fish populations, as well as the marine biodiversity23,27. Following the approach developed in ref. 5, extreme events in temperature (oxygen concentration) are defined as a prolonged instance with anomalously warm (low-oxygen) water (See “Methods” for the detailed definitions). We expect our findings to further deepen the understanding of the synergy of extreme events in the oceans.

Results
Our results suggest that regions with relatively high frequency of heatwaves tend to suffer a highly frequent occurrence of low-oxygen extreme events (Fig. 1). The ridgeline plots present the observed temperature and oxygen concentration anomaly since the 1960s, from which we can see significantly negative correlations between the two variables in most ocean regions. The correlations between temperature and oxygen concentration anomalies become less negative yet remain evident after data detrending, as detailed in Supplementary Note 1. Outputs from a climate model provide a more detailed and clearer evidence: the joint distributions of the two type of extremes show apparent overlaps between regions with marine heatwave and with low-oxygen extreme events, which indicates the spatial patterns of annual number of event days are highly coherent across the two types of extreme events (Fig. 1a, b). This relationship could also be found in the global patterns of individual extreme events (Supplementary Fig. 1): the pattern of low-oxygen extreme event frequency and duration broadly resembles that of marine heatwave events, with pattern correlation coefficients over 0.8 (p < 0.01). More importantly, both of the observations and model simulations reveal the impact of climate change on these extreme events is becoming progressively evident at the global scale, with an approximate increase in 0.07–0.09 months per year (Supplementary Fig. 2). In the late 20th century, only 10% of the world’s oceans are subject to ‘occasional’ occurrence of both heatwave and low-oxygen extreme events, but by the early 21st century they have expanded to 24.5%; conversely, over half of the global ocean (~69.3%) is under ‘rare’ occurrence of the two extreme conditions during period 1985–1999, while it drops sharply to about one quarter after year 2000 (Fig. 1c). The apparent connections between the two types of extreme events can be found in both duration and frequency, visualized through purple-colored regions in Supplementary Fig. 3. Most proportion experiences an increase in both two types of the extreme hazards rather than increase in one of the single hazard type, with hotspots of largest frequency increases occurring in subtropical regions of the North Pacific Ocean and high-latitude North Atlantic Ocean (Fig. 1a, b).

All of the above results illustrate the urgency and importance of studying the synergy of the two extremes. We thereby refer to the simultaneous occurrence of extreme high temperature and low oxygen conditions in space and time as a “compound low-oxygen extreme and heatwave” event (CLH event, see “Methods”). For the daily-mean model outputs, the duration for each CLH event is calculated in each grid cell and then used as a basis to figure out the annual statistics such as the frequency of the events, annual mean duration and total number of CLH days per year. The analysis of observation-based data is limited to basin-scale and the monthly time step due to a relatively insufficient spatial-temporal coverage. Furthermore, the likelihood multiplication factor (LMF, see “Methods”) is defined here to quantify the likelihood that a day (or month) is under CLH conditions relative to the likelihood assuming heatwaves and low-oxygen extreme events would occur independently from each other31–33.

The features of the LMF reveal statistical dependencies between marine heatwaves and low-oxygen extremes. The observational data indicates a globally averaged LMF of approximately 4.0. This implies a strong positive correlation between the two type of extreme events and, likewise, that the co-occurrence of the two extremes (i.e., CLH event) would be four times more frequent than the situation expected if heatwave and low-oxygen extreme were statistically independent. For the basin scale, the observed LMFs are larger than 1 in almost all of the ocean regions except the Eastern Pacific (colored dots in Fig. 1d). Model simulations are in good agreement with the observations, showing that over 90% of the ocean experience more CLH events (i.e., LMF > 1, contour lines in Fig. 1d), with a global average LMF value of approximately 5.1. Furthermore, the simulations reveal distinct LMFs among the global ocean due to differences in their intrinsic driving mechanisms. The largest LMF (> 5) occurs in the subtropical regions such as mid-latitude Pacific and Atlantic Ocean, indicating strong correlations between the two extremes in these regions, while the LMF < 1 in the equatorial East Pacific and Indian ocean in which relatively less probability of CLH events occurrences exist. We further perform the LMF analysis after detrending the data, serving as a sensitivity analysis within our study. The details can be found in Supplementary Note 1, illustrating the impact of temperature and oxygen concentration trends on the LMF of CLH events and demonstrating the robustness of our results.

The evolutions of CLH event reveal substantial modifications as a consequence of climate change (Fig. 2b). The annually observed number of months with CLH events has significantly risen in the majority of the oceans (Supplementary Fig. 4). Comparisons between two 15-year periods (1985–1999 and 2000–2014) in the climate model also show a considerable rise in annual CLH days, with increases in more than 70% of the world’s ocean and a global average increase exceeding 20 days, despite large spatial variations exist (Fig. 2c). The largest enhancement occurs in the mid-latitude ocean basins, with an increase of up to 100 days. Decreases in annual CLH days are found in parts of the eastern tropical Pacific and high-latitude Southern Ocean. As a global average, the annual CLH days have increased significantly in the past few decades, with a linear trend of +0.96 days per year (p < 0.01; Fig. 2b, purple line). By the 2010s, the global ocean increases by approximately 25 additional CLH days per year compared with a baseline level of 14 days in the 1980s. This increase in annual CLH days indicates either an increase in the frequency of CLH, a prolongation of event duration, or an increase in frequency in conjunction with a prolongation of duration. Supplementary Fig. 5 shows the characteristics of CLH frequency and annual mean duration, from which we can see the distinct contributions of the two variables to the increases of annual CLH days in different regions. The detailed discussions about the contribution of frequency and duration changes to the increase of annual CLH days can be found in Supplementary Note 2.

It should be specifically noted that, alongside the increasing number of event days, the low-oxygen extreme events are becoming increasingly associated with heatwaves over the past few decades: during the period 1970–1984, about 41.0% of the low-oxygen extreme days are accompanied with heatwave globally, while it increases to 45.5% in 1985–1999 and eventually exceeds over 50% in 2000–2014 (Venn diagram in Fig. 2b). In other words, there is one out of every two low-oxygen extreme days that simultaneously suffer extremely high temperatures in the early 21st century. At the regional scale, the strengthening of coupling can also be found in most oceanic basins except for the East Pacific Ocean, although slight difference exists in the magnitude (percentages in Fig. 2c). The joint occurrences (i.e., CLH events) have become the dominant concern in the evolution of low-oxygen extreme events, reflecting a substantial increase in its synchronization with heatwaves.

The increasing co-occurrence of heatwave and low-oxygen extreme could have far-reaching implications for human society. Here we present the characteristics of CLH events in four representative fishing waters (i.e., The North Sea, Hokkaido fishing zone, East China Sea and Newfoundland fishing zone, Fig. 3a–i). Comparisons between two periods (1970–1984 and 2000–2014) reveal that these regions suffer a significantly longer persistence of both extreme high-temperature and low-oxygen conditions under climate change (Fig. 3a, c, f, h). There has also been an apparent enhancement of extreme event intensity in these fishing zones over the last few decades. The associated temperature and oxygen concentration anomalies show that CLH events after 2000 bring extreme conditions of much higher temperatures (+13%) and lower oxygen concentrations (−11%) than those before 2000 (Fig. 3b, d, g, i). It is also noteworthy that the response of CLH events to climate change is more sensitive in regions with abundant fish biomass than in other marine areas. In combination with global ocean biomass map34,35, we found that the most substantial average increase in annual CLH days occurs in ocean regions with high fish biomass (Fig. 3j). In regions where fish biomass exceeds 50 g m−2, there is an upward trend of CLH days about 1.47 days per year. The rise is approximately 50% greater than the global average increase. Specifically, the increase in annual CLH days in the North Sea, Hokkaido, East China Sea, and Newfoundland fishing zones is +3.64, +2.95, +2.66, and +3.34 days per year, respectively, all of which are more than 2.5 times the global average (approximately +0.96 days per year). Only 8% of the ocean where fish biomass below 5 g m−2 suffers an increase of CLH days larger than the global average level. In contrast, this percentage rises to 66% in the ocean regions where fish biomass exceeds 200 g m−2. Fish biomass declines more rapidly with increasing trends of CLH days under climate change (Supplementary Fig. 6). This greater increase of CLH events in high-biomass regions implies more severe damage to marine life and ecosystem, which may eventually influence human communities that rely heavily on fisheries, especially for those coastal regions in developing countries. It should be noted that low-resolution climate models such as CESM are recognized for exhibiting biases in coastal regions due to their limited ability to capture small-scale processes such as boundary currents, coastal dynamics, and ocean eddy fluxes36. Previous study has revealed an underestimation of marine heatwaves in CESM within coastal regions throughout the historical period37. In other words, the actual conditions in coastal regions might be more severe than anticipated in our manuscript. This underscores the pressing need for high-resolution models operating at daily time-scales to explore extreme events in biomass-rich coastal oceans.

A systematic assessment of the drivers behind these changes is necessary and urgently needed. For this purpose, here we employ the ensemble empirical mode decomposition (EEMD), an adaptive one-dimensional time series analysis method that separates scales naturally without any prior subjective criterion38 (see “Methods” for detailed descriptions). As demonstrated in previous studies39–41, the EEMD method can effectively separate anthropogenic forcing signal from natural internal variability in the time series. For example, based on this method, we can split the time series of global averaged seawater temperature as well as oxygen concentration into relatively short-lived and transient oscillation component associated with natural internal variability and long-term warming signal primarily related to anthropogenic forcing, as depicted in Supplementary Fig. 7. For multi-dimensional spatial–temporal temperature and oxygen concentration data, the EEMD method is applied into each grid cell to obtain the long-term signal as well as the oscillation component for these variables. Following ref. 39, We then piece together the oscillation components from all grid cells to create a gridded dataset which constructs a temporal evolution of the spatially coherent structure of this component. Eventually, we can obtain the spatiotemporal evolution of extreme events related to natural internal variability by conducting extreme event detection on this dataset. Similar methodologies are employed on the long-term signal to investigate the impacts of anthropogenic forcing. By using this EEMD-based method, we can thereby obtain the characteristics of extreme events under influence of natural internal variability and long-term signals.

Figure 4a illustrates the annual changes in CLH days associated with the long-term warming and oxygen decline caused by anthropogenic forcing. It can be seen that the vast majority of the world’s oceans are facing an increase in the number of CLH days since the 1980s, with some regions experiencing an increase of more than 50 days (Fig. 4a). This spatial pattern is directly shaped by the concurrent trends of warming and decrease of oxygen concentration. For example, certain regions in the North Pacific exhibits the highest values, indicating that these oceans experience the most pronounced concurrence of temperature increase and oxygen concentration decline.

The CLH events related to relatively short-lived and transient oscillation component can be caused by the interaction of many local and remote processes and phenomena acting across a large range of temporal and spatial scales. The underlying mechanisms are highly intricate. The strength of many of local processes (e.g., heat fluxes and vertical mixing) are a function of the overlying atmospheric synoptic conditions. These conditions fluctuate in response to shifts in the general circulation and internal climate variability, such as El Niño–Southern Oscillation (ENSO) and Indian Ocean Dipole (IOD), across the globe. Therefore, the contribution of natural internal variability significantly varies by location. It has a positive impact on annual CLH days in widespread North Pacific and North Atlantic oceans during period 2000–2014, while negative impact is found in regions such as Equatorial Pacific and Southern Ocean (Fig. 4b). Based on method developed by ref. 42, we furtherly explored the occurrence of extreme events in relationship with various modes of internal climate variability which have characteristic time scales ranging from sub-seasonal (e.g., the Southern Annular Mode, SAM) to multidecadal (e.g., the PDO), as these climate modes act to modulate the local conditions - either from local sources or remote sources via teleconnection mechanisms. Their relationships are quite complex, but we can still find clear patterns consistent with the known temperature or pressure patterns (Fig. 4c; Supplementary Figs. 8, 9). For example, ENSO seems to be the main modulator of CLH in the East Equatorial Pacific as well as Northwest Indian Ocean; the SAM dominates the occurrence of CLH event in vast ocean south of 40 S; the IOD, represented by Dipole Mode Index (DMI), strongly relate to CLH events around Indonesia. The region of North Atlantic in 30–50°N is very interesting that this region is influenced by various climate modes, including NAO (Supplementary Fig. 10). However, the influence of NAO is relatively minor compared to other modes; hence it is not drawn in the figures.

Our results reveal that both anthropogenic forcing and natural internal variability are responsible for the modifications of CLH events: the former is the main contributor to the lengthening of annual CLH days at the global scale, while the latter plays an important role in spatial distributions of this change. For instance, mid-latitude oceans, including the North Pacific and North Atlantic, experience substantial influence from the long-term warming signals driven by anthropogenic forcing, as well as variability modulated by climate modes such as PDO (Fig. 4). The collective impact of these factors has contributed significantly to a noteworthy increase in CLH days within mid-latitude oceans. It’s worth noting that the signals separated by the EEMD method are generally considered to be caused by anthropogenic forcing and internal climate variability. We acknowledge that our method does not identify specific signals such as aerosols and volcanic activities. Further attribution analysis is needed to disentangle the contributions of various factors in the future.

A range of mechanism can lead to coherent synergy of anomalously warm and low-oxygen waters. Temperature variations may result from a confluence of processes, including air-sea-air heat exchanges, current and eddy-driven advection, horizontal and vertical mixing, as well as the entrainment of water into the surface layer19. For the oxygen concentration, changes in oxygen solubility due to sea surface temperature variations plays a crucial role in modifying these awater’s capacity to store oxygen. Concurrently, the air-sea oxygen flux, along with the horizontal and vertical transport of dissolved oxygen as well as local processes of oxygen production and consumption, collectively shape surface dissolved oxygen concentrations. These processes mentioned above are subject to modulation by local and remote factors operating over a broad spectrum of spatial and temporal scales. The strength of many of these processes are a function of the overlying atmospheric synoptic conditions. These conditions fluctuate in response to shifts in the general circulation and climate variability, such as ENSO, PDO, and IOD, as well as the human-induced long-term changes across the globe42. Under the combined dominations of natural internal variability and anthropogenic climate change, These modifications lead to the accumulation of heat and eventually outbreak of heatwave19,43. This extreme high temperature simultaneously reduces the solubility of dissolved oxygen, which directly cause oxygen escape from ocean10,44. For further validation, we investigated the changes in oxygen solubility during extreme events (Supplementary Fig. 11a). The oxygen solubility in seawater during CLH events exhibits an apparent decrease compared with the average state, indicating a reduced capacity of seawater to store oxygen. Interestingly, enhanced stratifications, indicated by mixed layer shallowing, occurs in the ocean during the events (Supplementary Fig. 11b). This intensified ocean stratification can hinder the ventilation of oxygen13, suppress the ocean’s oxygen supply, which might eventually affect oxygen in the ocean interior.

Methods
Dissolved oxygen and temperature observations in the ocean
We use in situ observations from World Ocean Database (WOD) 2018, the world’s largest collections of the uniformly-formatted and quality-controlled ocean profile data47. Approximate 1.3 million profiles flagged good are selected from the database, which contain both temperature and dissolved oxygen data since 1960 (Supplementary Fig. 12). Given our primary emphasis on extreme events in the uppermost ocean layer, only profile data near the ocean surface (0–10 m) is utilized in this study. We referred to Schmidtko’s paper4 to complete a series of technical steps such as data filtering and preprocessing (quality control, duplicate removal, etc.), except that the ocean is divided into nine basins for sake of comparisons and analysis (the definition of ocean basins can be found in Supplementary Fig. 13). We first divided the ocean into grids of 2.5 degrees by 2.5° and mapped each profile data to the grid point closest to the recorded latitude and longitude at the time of observation. Then, we followed the defined basin boundaries, performed a weighted average of all eligible data based on their grid coordinates, and eventually obtained the time series for each basin. Note that this observational data is limited to monthly-mean resolution due to insufficient temporal coverage. We acknowledge that the averaging in the mapping is inevitably affected by missing values. However, limited by the spatiotemporal coverage of WOD data, this is the compromise approach we adopted after comprehensive consideration. The World Ocean Atlas (WOA) 2018, which is based on WOD18 and provides climate mean map of the temperature and dissolved oxygen, is used in this study as the reference to evaluate the model simulations.

Two additional observation-based datasets, namely the Optimum Interpolation Sea Surface Temperature (OISST)48 dataset and the Gridded Ocean Biogeochemistry from Artificial Intelligence-Oxygen (GOBAI-O2)49 dataset, are also employed in this study. The OISST dataset is derived from remotely sensed sea surface temperature by the advanced very high-resolution radiometer infrared satellite data and in situ measurements. It maintains a global spatial resolution of 0.25° on a daily scale since 1981. The recently published GOBAI-O2 dataset provides ocean dissolved oxygen concentrations derived through machine learning algorithms trained on oxygen concentration observations from Argo-float-mounted sensors and discrete measurements from ship-based surveys. It covers the period from 2004 to 2022 with a monthly 1° × 1° resolution. For our study, both the OISST and GOBAI-O2 datasets are processed to a consistent spatial and temporal resolution and specific overlapping time periods when both datasets are available.

Descriptions of the climate model
The global analysis of extreme events was mainly conducted with outputs derived from the fully-coupled Community Earth System Model (CESM, version 2.1.3) developed at National Center for Atmospheric Research (NCAR)30. The CESM consists of ocean, atmosphere, land, sea-ice, land-ice, river, and wave models, which could provide state-of-the-art computer simulations of the Earth’s climate states. For the ocean component, the description of physical process is based on the Parallel Ocean Program version 2 (POP2) and the ocean biogeochemistry is simulated by Marine Biogeochemistry Library (MARBL). The ocean has a normal 1° × 1° horizontal resolution (gradually increasing toward the Equator), and 60 vertical levels. In this study, Our CESM experiment is the same as the CMIP6 historical experiment50. As a result, the forcings are based on observations from 1850 to 2014 and include greenhouse gas emissions, land-use conditions, solar forcing, and so on. The forcing data, in particular, were officially provided by the CESM group, as introduced by Danabasoglu et al.30 in Section 5. Only the last 55 years of the historical experiment (i.e., 1960–2014, period when the intensity of human activity increases dramatically) is used for validations and analysis.

Definition of extreme and compound extreme events
A seasonally varying climatological threshold is used to define the extreme event, that is, the temperature above 90th-percentile threshold for marine heatwave (extremely high-temperature event)15,51 and the oxygen concentration below 10th-percentile threshold for low-oxygen extreme5,21. In this study, the percentile threshold is calculated over the full range of the time series studied in our manuscript, that is, from year 1960 to 2014, and determined from the shortest available time-step, i.e., monthly for observations and daily for the model simulations. For each grid cell or ocean basin, we separately calculate the threshold for each calendar day (model outputs) or month (observation data) of a year. To ensure the continuity of extreme events, for the daily time series, a duration of at least 5 days with temperature above the threshold (oxygen concentration below the threshold) is imposed in the identification of an extreme event, and two events separated by an interval of two or fewer days are considered as one event.

The compound low-oxygen extreme and heatwave (CLH) events are defined as the days or months when extremely high-temperature and low-oxygen conditions occur simultaneously (that is, both temperatures above 90th-percentile and oxygen concentration below 10th-percentile at the same time and location). Similar to marine heatwave and low-oxygen extreme events, we also impose a criterion on the minimum duration of this compound event, with the caveat that it is slightly shorter compared with the individual extreme event, specifically 2 days. The annual statistics for CLH events, such as event frequency, annual mean duration, and the total number of extreme days per year, are calculated in the same manner as marine heatwaves.

Evaluation of model simulations against in situ observation
Comparisons with WOD18 and WOA18. We firstly evaluate the simulated oceanic dissolved oxygen and temperature in comparison with the observational data from WOD18. As shown in Supplementary Figs. 14, 15, the model is able to capture the observed temperature and dissolved oxygen patterns. For the climate mean state, the model simulates well the distribution of temperature and dissolved oxygen concentration among the ocean basins. It should be noted that, similar with results from other CMIP-type models36, bias exists in the our CESM outputs. A Lower temperature, about −0.42 °C on average, is simulated in some high-latitude regions compared with observations, while the opposite is found in part of central-eastern Pacific and Atlantic, approximate +0.39 °C on average. The simulated dissolved oxygen exists bias mainly in Extra-equatorial ocean, where the model overestimates oxygen concentration with an average of approximate +4.69 mmol m−3. Nevertheless, the simulated distribution is in good agreement with it in observations from WOA2018, with pattern correlation coefficient over 0.90, significant at 0.01 level. More importantly, the model captures very well the temporal evolution of the extreme and compound extreme events under climate change (Supplementary Fig. 2). The marine heatwave is simulated to increase by 0.09 ± 0.01 month per year since the 1960s, which corresponds well with the observed trend of 0.07 ± 0.01 month per year. Besides, both the model and observation reveal an increase in low-oxygen extreme events, approximate 0.09 ± 0.01 month per year over 1960–2014. Accurate simulations of the heatwave and low-oxygen extreme events eventually lead to a very good agreement between simulated CLH event and the observation-based characteristics.

Comparisons with OISST and GOBAI-O2. We also assess the model simulations against the results from OISST and GOBAI-O2 datasets. Supplementary Figs. 16, 17 include climatology comparisons between the simulation and observation. This analysis spans the timeframe when data are concurrently available in both observation datasets and model simulations. The bias and root mean square error (RMSE) metrics for sea temperature reveals that, overall, the simulation performs well, with biases typically less than 1 °C. Noteworthy exceptions, such as Kuroshio, North Atlantic coast, equatorial South Africa, Peru coast, California coast, have been identified and discussed in earlier studies. The distribution of RMSE is consistent with that of bias. The larger bias corresponds the large value of RMSE. Therefore, the simulation has a good performance to capture the variation of sea surface temperature. Supplementary Fig. 17c illustrates that the bias for oxygen concentration is generally less than 15 mmol m−3, representing approximately 10% relative to the climatology. Despite larger biases observed in specific regions, such as the Indian Ocean and high latitudes of the Atlantic and Pacific Oceans, we think that these discrepancies fall within an acceptable range for investigating oxygen variation. The distribution of RMSE between simulation and observational aligns with that of bias, with the majority of values being less than 20 mmol m−3. We consider such bias acceptable for investigating the variation in oxygen levels.

The trend distribution of Seawater temperature and oxygen concentration is illustrated in Supplementary Figs. 18, 19. In general, the simulated temperature trend aligns well with the observations, although the trend of seawater temperature in the North Pacific is overestimated, and in the South Pacific, it is underestimated. The bias of trend in the Atlantic is much less than that in Pacific. Regarding oxygen concentration, the simulations indicate an overall decline in most ocean regions, which is consistent with the observations, despite that the simulation overestimated the decreasing trend in the Pacific and underestimated the increasing trend in the Atlantic, and in the southern hemisphere, the increasing trend in dissolved oxygen concentration contrasts with the observed decreasing trend. The temporal correlation coefficients between the observed and simulated seawater temperature, as well as dissolved oxygen concentration, were investigated in Supplementary Fig. 20. The results suggest that the correlation coefficient exceeds 0.9 in most regions, except for certain oceanic areas near the equator. This finding highlights a noteworthy coherence between the observed time series of seawater temperature and oxygen concentration and those derived from simulations.

We explore the spatial distribution characteristics of CLH events using observation-based data and compare them with model-derived simulations, as depicted in Supplementary Fig. 4. Notably, the time resolution of GOBAI-O2 is monthly, thus the unit for observation-based annual CLH days is denoted as “months per year”. For the purpose of comparison in Supplementary Fig. 4, we convert this unit to “days per year” by assuming a factor of 30 days per month. Due to limitation of temporal coverage of GOBAI-O2 dataset, the trend of annual CLH days for observation is calculated for period 2004–2022. The simulated climatology of annual CLH days exhibits remarkable similarity to that derived from observations, with a pattern correlation of 0.84. Regions with the most pronounced annual CLH days are situated in the North Pacific, mid-latitude North Atlantic, and certain areas in the South Indian Ocean. Despite existing differences, both simulations and observations indicate a global increase in annual CLH days. A significant upward trend is found in the North Pacific and the vast Atlantic Ocean in both simulations and observations. Increases in annual CLH days are also evident in regions of the Indian Ocean and Southern Ocean, although the rise in simulations is generally less pronounced than in observations.

The comparisons between the high temporal resolution model simulations and observational data also suggest that the feature of extreme events is relatively insensitive to the temporal resolution. The great consistence between model and observation gives us confidence in using the climate model (CESM) to analyze the spatial and temporal evolution of ocean biogeochemical extremes and heatwaves at the global scale.

Classification of the extreme events
To present a clear description of extreme events and their joint distribution, for each type of extreme event, we firstly classified the number of event days into four levels, that is, so-called Rare (R), Occasional (O), Frequent (F), Common (C) in our study. The thresholds dividing these four levels are the 40th, 70th, and 90th percentiles of annual event days across the oceans (Supplementary Fig. 21). After classification, the evolution of extreme events and their joint distribution are straightforward in Fig. 1a–c.

Likelihood multiplication factor
The LMF is used to illustrate the impact of the possible correlation between hazard pair (that is, marine heatwave and low-oxygen extreme) on their joint occurrence probability32, and likewise, estimate in which regions the compound events are more likely to occur33. The LMF in this study is given by:

LMF = P(cooccurence of heatwave and low oxygen extreme) / (P(heatwave) * P(low oxygen extreme)) (1)

the ratio of actually observed probability of joint occurrence (i.e., CLH event) and probability assuming the two types of extreme events complete independence.

Depending on the relationship between the hazards pairs, the LMF can range from 0 to an upper limit defined by percentile thresholds (it is 10 in our manuscript). An IMF higher than 1 indicates compound event occurs more often than by chance, and likewise, that there is a positive dependence between heatwave and low-oxygen extremes, with the value increasing with the strength of correlation. In contrast, LMF less than 1 suggests that the two types of extreme are negatively correlated, that is, a reduced likelihood of compound event. The LMF is calculated for each grid cell (model outputs) and each ocean basin (model outputs & observations) in this study.

Ensemble empirical mode decomposition method
The EEMD is an adaptive one-dimensional time series analysis method that separates scales naturally without any prior subjective criterion52. This noise-assisted technique could perform operations that partition a series into different ‘modes’:

X(t) = sum(IMF_i(t)) + r(t) (2)

where IMF_i is the ith Intrinsic Mode Function (IMF) which are oscillatory components with distinct frequencies and amplitudes and r(t) is the residual part of the data after all the IMFs have been extracted.

In this study, the EEMD is used to separate signals of anthropogenic forcing and natural internal variability in the evolution of extreme events. Time series of sea temperature and oxygen concentration are decomposed via the EEMD method for each grid cell. We set the amplitude of added noise to 0.2 times the standard deviation of the original data. The ensemble number is 400 and the total number of IMFs is 6. Following ref. 40, the sum of IMF 1–5 is treated as the relatively short-lived and transient oscillation components related with natural internal variability and IMF 6 is the long-term signal for anthropogenic forcing.

Using the EEMD method, the original time series is decomposed into two components: the anthropogenic forcing signal and the internal natural variability signal. By performing extreme event detection on the anthropogenic forcing signal, we can identify extreme events related to anthropogenic factors. Similarly, extreme events related to natural internal variability are identified by conducting extreme event detection on the latter.

An important clarification can be found in Supplementary Note 3 to perceive the underlying consistency between the results from EEMD-based method in our study and the shifting baseline approach suggested by Amaya et al.53. The principles underlying the EEMD-based method resonate with the fundamental concept articulated by Amaya et al. (2023), which highlights the need for separation of influence of long-term signal and relatively short-lived and transient changes on extreme events.

Attributions to large-scale modes of climate variability
An analytical framework established by Holbrook et al.42 is used in this study to identify regions where there is a statistically significant relationship between the occurrence of CLH events and the internal natural variability represented by large-scale climate modes. This framework provides a systematic approach for evaluating the impact of climate variability on extreme events, and can be applied to a wide range of data and models.

Nine commonly used climate indices are considered in this study, including Niño-3.4 index, EMI (ENSO Modoki Index), DMI, NAO (North Atlantic Oscillation) index, PDO (Pacific Decadal Oscillation) index, TPI (Tripole index), Atlantic Niño index, SAM (Southern Annular Mode) index, NPGO (North Pacific Gyre Oscillation) index (see Supplementary Note 4 for detailed descriptions). For each grid, we calculate the frequency of compound events during positive and negative phases of these key climate modes and identify which climate modes are significantly related to occurrence of CLH event following the method by Holbrook et al42. In our manuscript, this framework is done by summing the number of CLH days at each location and determining whether the climate index is in a positive or negative phase. Monte Carlo simulations are conducted to determine whether the number of days counted were greater or less than what might be expected by chance. For a specific climate index, we generate a synthetic time series with autocorrelation characteristics akin to the index, utilizing a fourth-order autoregressive model. Subsequently, we recalculate the number of CLH days during positive and negative phases of the synthetic index. This process is repeated 10,000 times to produce a frequency distribution of the expected number of days. The 5th and 95th percentiles are used to form confidence intervals for the given mode and region. The climate modes (and phases) which have the greatest significant impact on occurrence of CLH event are eventually shown in Fig. 4c. By using this framework, we aim to gain a better understanding of the underlying causes of CLH events linked to internal climate variability.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Review of the “The ocean losing its breath under heat waves” by Changyu Li, Jianping Huang, Xiaoyue, Lei Ding, Yongli He\nTitle: I like it.\nAbstract: please make clear that this is only considering the surface ocean (1-10m) from the outset.\nGeneral Comments:\nI am quite happy with the approach of the EEMD. The long term trend is not necessarily human induced change, because that would require all of the forcings to be separated out (ie CO2 and aerosols and also the natural variations like volcanic eruptions and then internal oscillations. You don't actually attribute these changes to changing forcings (in the way that the attribution literature and IPPC do). I agree this is beyond the scope of this work. I would simply say the trend found at each grid point is a long term climate signal likely to be attributable to changing forcing , and leave it to that.\nOverall I think this is a nice piece of work. It connects oxygen waves and heatwaves and their compounding in surface waters. It does it with as much observations as is available and also with simulations. Overall it is in my view a compelling study (in spite of my comments below). Oxygen decline in the oceans, in extremes and also as trends in the surface ocean has a lot to say about changes in surface biology and stratification. This issue needs to be brought forward to the wider community.\nMinor comments:\nLine 29: to suffer highly frequent occurrence\nLine 92: .. with an approximate increase in duration of 0.07 …\nLine 146: why the largest increase in mid-latitude ocean basins…? Explain in “discussion section”\nLine 150: … the global ocean increases by approximately 25 additional …\nLine 227. Delete the phrase “induced by human activities.” for the reasons explained above.\nLines 249- 258: I agree these patterns are very interesting. Is this relationship to the modes based on the time series of long term trends (described earlier) or from the time varying trends (described earlier), and the indices used for the different modes, are they derived from the model simulations or somewhere else?\nLines 266-268: I like this conclusion. The lengthening is also found in other studies of marine heat waves for example.\nLines 269-296: This is the explanation of the processes. Now that I understand this is for surface waters (0-10m) layer it is easier to understand the outcome of the results. Although this text addresses my original review, it doesn't actually address it in any quantitative way. In the first review I was curious to know if the extreme heat waves were affecting the ocean deoxygenation by the amount expected from the solubility decrease, or is it larger in some way and thus revealing something about productivity of the region. The LMF factor shows that events co-exits and compound but doesnt show it relative to the solubility change in oxygen driven by a warm event. Did the onset of heat wave alter the stratification and thus make it more difficult for oxygen to be ventilated in these surface waters. The relationship to these two factors (solubility and stratification) would be interesting in terms of explaining the compounding of heat waves with low oxygen events, and more useful that the current text that talks about the mechanisms in these lines..\nSpecific Comments: Every figure must be clear on the depth level that is being analsyed or compared. Every figure in the responses to the authors never mentions the depth layer that is being shown. Depth is never mentioned in the responses text. There is only one place in the main text where the depth is revealed, and that is well into the methods section, and yet is should appear in every figure caption. Please introduce the surface depth range and define the layer of interest very early on in the text, before line 64, and in the abstract (page 1). Every figure caption should say surface (0-10 m). Much of the oceanography community think about the full water column, whereas this paper is only about the 0-10m surface layer of the ocean. Makes a huge difference in communication.\nFigure R8. Once again this figure should say surface (0-10m) layer or similar. Is this the trend from the beginning of the record (1960 to 2014) say? The panels of this figure are showing that the correlation with the time trend in the data.\nI was and remain confused on what is presented in the figures and what is a result from a simulation, and what is a result from an observation? While it could easily be stated in the text, it helps enormously to know from the caption whether a particular panel is composed of observation an or simulations. See the questions about each figure.\nFigure 1: This paper uses both models and observations. In this figure I cannot tell which panel is an observation? Which panel is a simulation. In panel d, I think this is simulation, but I cant see dots that are observations? I can see histograms/ridgelines connected to dots on the map. Are histograms/ridgelines composed entirely of observations of temperature and oxygen in the regions whose boundaries are explained in the response to reviewers? Are they simulations?\nFigure 2: This panel is clear that each panel is simulations. Results entirely from the ensemble of 5 models.\nFigure 3: This is a complex figure. Are a-d observations or are the simulations? Or is there a mixture of simulations and observations in these panels? I can't tell from the Figure caption. Same question for panels f-i as for a-d. Panel e: is this also simulated fish biomass (eg from FISHMIP) or is it an observation? Panel J: not sure how to interpret panel J, does it imply higher biomass is associated with longer CLH days? While it shows big CLH are associated with high biomass, it doesnot show that areas of high biomass are more affected by high CLH days than low biomass, and isnt this the more important result but shown in Fig S6? Compare with text in lines 183:185Lines 192-193. Panel k: Why are there 6 regions (pie charts) when there are only 4 biomass locations shown in panel e? Is Panel k observations or simulations?\nFigure 4: Panel a: “anthropogenic component” is based on simulations (but through the EEMD estimation. Panel b: is also based (and estimated from) on simulations. It still needs stating. Units should be added to the colorbar. I like panel C. I am assuming this is from simulations (analysed 5 times from the 5 ensemble simulations) and the 9 modes climate are each diagnosed from each model simulation using the model indices? Is that right? I like the map in C and the bar chart. Is this map (and bar charts) estimated from the panel b natural variability or from the time series of both panel a and panel b. What is the period for this map, I find it a bit odd to talk about increased occurrence of CLH and connecting that directly to the dominant modes. Does this include anthropogenic forcing or not? And for what time period is the enhancement of occurrence estimated over, and for all of the ensemble runs?", "reviewer_2": "The authors have carefully revised the manuscript in response to the reviewers comments and have added interesting and helpful analyses. In particular, I think the additional supporting information analyses, such as comparing EEMD results to those based on a shifting baseline, are very interesting. Overall, the methods underlying the results are now much clearer to me. Despite these improvements, I still have remarks and questions that should be addressed before publication.\n##########\nMajor comments:\n##########\nI am somewhat missing a more careful embedding of the present study in the published research on low oxygen extremes and their interaction with marine heatwaves in the introduction section. I think this point is very important. Otherwise the contribution of this paper, what it adds to existing knowledge on LOEs and CLH events, can not be properly put into context. See more specific comments below.\nI am a bit confused about the identification of correlations and definition of the extremes from the WOD data (e.g., Supplementary Figure S2) across the ocean basins, or more general how the WOD data is handled. While you nicely explain how the extremes are defined from the monthly data in the 9 basins, it is still illusive to me how the data is aggregated over the 9 basins in the first place. In the methods, you refer to the Schmidtko study - however it is not explained yet what is done exactly. Do you just average all measurements in a month, say November 2008, over a whole basin, to eventually get one monthly time series for each basin? If this is done, then it is questionable how this procedure influences correlations and compound event occurrence. For example it may be that in a November, only measurements from the gulf of Alaska enter the time series of the North Pacific, with much different mean temperature and oxygen compared to the subtropical North Pacific. My thoughts here may not apply, but a clarification is needed.\nLow oxygen extremes are often looked at in the ocean interior. I think the authors need to establish why low oxygen extremes at the surface are relevant to make strong claims such as in L71-L74. Can you give a reference for this claim? Likewise, the final sentence of the manuscript, 'Under their combined effects, marine heatwave amplifies the peak signal of deoxygenation, highlighting urgent needs for adaptation responses planning for future CLH risks.', is too strong if the importance of low surface oxygen events can not be established.\n##########\nSpecific comments:\n##########\nL34-35: \"and state-of-the-art climate model simulations\": given that you, as far as I see, only analyze one climate model simulation, this should be \"and a state-of-the-art climate model simulation\"\nL39: \"Our work shows a widespread emergence of low-oxygen extreme events in response to increasing marine heatwaves\" -> as you write in the response to reviewers, and have already corrected earlier in the abstract, isn't it rather mainly an increase in low-oxygen extreme events due to increasing mean temperatures than due to an increase in marine heatwaves?\nL58: I would propose adding Hobday et al., 2016 and Frölicher et al 2018 to the references on marine heatwaves\nL58-59: There is some work on low oxygen events (as you also write before) and also to a lesser extent on CLH events, for example Köhn et al on oxygen extremes, this should be added here as an example of studying low oxygen extremes, and also Gruber et al. 2021, as well as Mogen et al 2022 (https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2021GL096938) who specifically looked at the heatwave impact on oxygen. Overall, I think the results of the current study needs to be more carefully embedded in existing literature. I would propose going through section 5.3 of this recent discussion paper https://bg.copernicus.org/preprints/bg-2023-182/ for further references for studies about low oxygen extremes. The knowledge gap needs to be stated more precisely.\nFigure 1: The corresponding correlation coefficients are presented in the top right corner with symbols * and ** representing the statistically significant at P<0.01 and P<0.1 level, respectively.\n-> * and ** were switched in this sentence, ** should refer to P<0.01 I think.\nL79-L82: The additional analysis on the correlation coefficient and LMF if the data is detrended (Text S1 and Figure S21) provides very important information. Later on, you discuss that the LMFs > 1 are weaker but still present when detrending the data (L135-139). I would propose to add also somewhere that the correlations of T and O2 anomalies (shown in Fig 1d) are less negative but still present after detrending the data.\nL83-86: Thank you for looking into the other metrics in response to my comments. I am not yet fully convinced by the reply in the R2R about why the event frequency is the best choice here. Given that you show a composite of the number of extreme days for single and overlapping events in the LMF map in Fig 1d, the most useful choice for Fig 1a,b would, in my opinion, be to also show number of event days here. The frequency of individual events is not only influenced by how frequently the variables are above the thresholds (i.e., the number of days), but also inversely controlled by how long individual events take on average, and thus somewhat harder to interpret.\nL93-95: Fig 1a,b is clearly interesting, but it does not directly provide proof about the claim you make here. You show that areas with increasing MHW number are often also areas with increasing LOE number. However, this does, strictly speaking, not imply that MHWs and LOEs co-occur more often, just that both events individually become more frequent in a region. You discuss co- occurrence starting L107ff and Fig. 2 - I would leave it out here where it is not yet shown.\nL121-122: How is the global LMF estimated from the WOD data - is this the area-weighted average of the colored dots in Fig 1? and what is the global mean LMF from the CESM model? this information could be provided in L129\nFigure 1: It would be nice if the maps could be all centered the same way. This would for example help identifying that the regions with high number of CLH days (Fig 2a) are also those with high LMF (Fig 1d).\nFigure 2: The caption should state which panels are observational (WOD) data, and which are from CESM.\nL144: what is considerable increase? this should be specific when also given a specifc ocean area (> 70%) for it.\nL161-164: I see that the CLH fraction of LOEs is increasing, however, the proportion of marine heatwaves that are also LOE does not seem to increase at the same time (47% at the beginning and at the end). This somewhat weakens the statement made here in my opinion.\nL180: the normalization in Fig 3b,d,g,i - are the two periods normalized independently from each other? this is relevant to interpret the percentage changes in the normalized anomalies.\nL185: \"reaching a peak of 2.53 days per year. This rise is approximately 2.6 times greater\" I think it would be more informative to compare mean CLH trends over high fish biomass regions with global mean treands, rather than comparing a peak to a the overall mean.\nL223f: The decomposition of the changes in CLH events shown in Fig 4a,b is very interesting. However, while it is clear to me how Fig 4b is calculated (remove the anthropogenic forcing signal from the time series using the EEMD approach), I do net yet understand how Fig 4a is calculated exactly. Some more explicit information would be good here. Maybe as a sentence appended to the methods section paragraph on EEMD?\nL247: 2000 or 1980? since referring to Fig 4b, it should be 2000 I think. Furthermore please add the time period analyzed for Fig 4c to the caption.\nL261f: a smaller thing: what is meant with \"not significant compared with other climate modes\"? Usually statistical significance would be assessed for each mode contribution individually I think\nL267: I would maybe also add that the former is a main contributor at the global scale (compared to the internal variability that is more important for the spatial patterns of changes)\nL275f: I agree that air-sea oxygen is most likely very important - however I think the explanation here may not be very clear yet. I assume the variations in surface oxygen content from air-sea oxygen exchange or not so much driven by variations in atmospheric oxygen but by changes in oxygen solubility, mainly due to changes in SST, causing a change in the capacity of the surface water to hold oxygen.\nL287: the suppression of oxygen supply through stratification applies to subsurface waters, but I don't think it is important for surface waters that are studied here. Please rephrase.\nL292-295: Do you show that the changes during the extremes are larger than those expected from decades of anthropogenic changes? This statement may be a bit strong and not hold generally - please check.\nL295f: You show some trend that LOEs tend to be more often also MHWs over time, but it does not seem like MHWs are more often also LOEs (my comment above), 'intensified synergy' is maybe thus a bit too coarse?\nSupporting information text S3: This addition is very useful. I just want to make one comment as just my perspective on it. A shifting baseline, in my understanding, does not necessarily need to be based on a shifting window for extreme event definition (that is discussed in this section). It can also rather be a shifting-mean baseline - already realized by simply detrending the data based on some linear or polynomial spline, or subtracting the ensemble mean to identify forced changes (e.g., Burger et al, 2022), some low-pass like filtering, or the EEMD approach you use. After all, I would consider the EEMD a (notably very interesting) option to identify the forced changes for applying a shifting-mean baseline.\n##########\nMinor comments:\n##########\nL79 and at other places: The authors write \"extreme event\", although it should be \"extreme events\" in the plural form, if I am not mistaken.\nL82: \"Outputs from climate model\" -> from A climate model\nL167: Typo \"although slightly difference exists in the magnitudes\"\nL189: \"(approximately +0.96 days per year)\" This information is already given before."}, "author_response": "We appreciate the two reviewers for the time and effort expended in reviewing the\nmanuscript and for the constructive comments, which are in blue italics below.\nResponse to Reviewer #1’s Comments\nTitle: I like it.\nResponse: Thank you for your kind words regarding the title of our manuscript. We\nappreciate your feedback and are encouraged by your positive response.\nAbstract: please make clear that this is only considering the surface ocean (1-10m)\nfrom the outset.\nResponse: We have revised the abstract to clarify this point, which you could find in\nLines 34-36 in our revised manuscript:\n“By utilizing in-situ observations and state-of-the-art climate model simulation, we\nprovide a global assessment of the relationship between the two types of extreme events\nin the surface ocean (0-10m).”\nGeneral Comments:\nI am quite happy with the approach of the EEMD. The long term trend is not\nnecessarily human induced change, because that would require all of the forcings to\nbe separated out (ie CO2 and aerosols and also the natural variations like volcanic\neruptions and then internal oscillations. You don't actually attribute these changes to\nchanging forcings (in the way that the attribution literature and IPPC do). I agree\nthis is beyond the scope of this work. I would simply say the trend found at each grid\npoint is a long term climate signal likely to be attributable to changing forcing , and\nleave it to that.\nResponse: Thank you for your insightful comments and acknowledgment of the\nmethods employed in our paper. Yes, as you pointed out, these changes include\ninfluences from other forcings. Since the primary aim of this paper is to explore the\nanthropogenic and natural impacts on extreme events, we did not engage in further\ndiscussion on other forcings, as you said it is beyond the scope of this work.\nNevertheless, some statements in our manuscript may be somewhat imprecise. We\napologize for any ambiguity in our wording. In the revised version, we have changed\n“long-term warming signal related to anthropogenic forcing” into “long-term warming\nsignal primarily related to anthropogenic forcing”. Additionally, we have added a\nparagraph in the discussion section to provide further clarification on this matter (Lines\n295-300 in the revised manuscript):\n“It's worth noting that the signals separated by the EEMD method are generally\nconsidered to be caused by long-term signal (mainly as a consequence of anthropogenic\n1\nforcing) and internal climate variability, though we acknowledge that some signals\n(such as aerosols and volcanic activities) are still mixed in. Further attribution analysis\nis needed to disentangle the contributions of various factors in the future.”\nOverall I think this is a nice piece of work. It connects oxygen waves and heatwaves\nand their compounding in surface waters. It does it with as much observations as is\navailable and also with simulations. Overall it is in my view a compelling study (in\nspite of my comments below). Oxygen decline in the oceans, in extremes and also as\ntrends in the surface ocean has a lot to say about changes in surface biology and\nstratification. This issue needs to be brought forward to the wider community.\nResponse: We are grateful for your thorough review and constructive feedback, which\nhave undoubtedly strengthened the quality and relevance of our study. Your recognition\nof the importance of oxygen decline in the oceans, resonates deeply with the objectives\nof our research. We totally agree that these issues warrant greater attention and\ndiscussion within the wider scientific community.\nMinor comments:\nLine 29: to suffer highly frequent occurrence\nResponse: Done.\nLine 92: .. with an approximate increase in duration of 0.07 …\nResponse: Thank you for your careful review. This sentence has been revised.\nLine 146: why the largest increase in mid-latitude ocean basins…? Explain in\n“discussion section”\nResponse: We think that the rapid warming trend of mid-latitude oceans\n(Supplementary Fig. S18), combined with the influence of various climate modes on\nthese regions (Fig. 4), has jointly led to the largest increase in CLH days in mid-latitude\nareas. The rapid increase in extreme events in mid-latitude oceans can also be found in\nprevious research (e.g. Burger et al., 2018, which investigated compound marine\nheatwaves and ocean acidity extremes). Following your suggestion, we have\nincorporated the discussion in Lines 291-295:\n“mid-latitude oceans, including the North Pacific and North Atlantic, experience\nsubstantial influence from the long-term warming signals driven by anthropogenic\nforcing, as well as variability modulated by climate modes such as PDO (Fig. 4). The\ncollective impact of these factors has contributed significantly to a noteworthy increase\nin CLH days within mid-latitude oceans.”\nBurger, F. A., Terhaar, J. & Frölicher, T. L. Compound marine heatwaves and ocean\nacidity extremes. Nat. Commun. 13, 4722 (2022).\nLine 150: … the global ocean increases by approximately 25 additional …\nResponse: Done.\nLine 227. Delete the phrase “induced by human activities.” for the reasons explained\nabove.\nResponse: Done.\nLines 249- 258: I agree these patterns are very interesting. Is this relationship to the\nmodes based on the time series of long term trends (described earlier) or from the\ntime varying trends (described earlier), and the indices used for the different modes,\nare they derived from the model simulations or somewhere else?\nResponse: This section analyzed how climate pattern such as PDO and ENSO modifies\nCLH event. It is based on the time series from the time varying trends. The indices used\nto characterize the different modes are derived from the model simulations. The method\nfor computing these indices is outlined in Supplementary Text S4.\nLines 266-268: I like this conclusion. The lengthening is also found in other studies\nof marine heat waves for example.\nResponse: Thank you for your favorable comment regarding our conclusion. We are\ngrateful for your acknowledgment of the consistency between our findings and those of\nother researchers in this field.\nLines 269-296: This is the explanation of the processes. Now that I understand this\nis for surface waters (0-10m) layer it is easier to understand the outcome of the results.\nAlthough this text addresses my original review, it doesn't actually address it in any\nquantitative way. In the first review I was curious to know if the extreme heat waves\nwere affecting the ocean deoxygenation by the amount expected from the solubility\ndecrease, or is it larger in some way and thus revealing something about productivity\nof the region. The LMF factor shows that events co-exits and compound but doesnt\nshow it relative to the solubility change in oxygen driven by a warm event. Did the\nonset of heat wave alter the stratification and thus make it more difficult for oxygen\nto be ventilated in these surface waters. The relationship to these two factors\n(solubility and stratification) would be interesting in terms of explaining the\n3\ncompounding of heat waves with low oxygen events, and more useful that the current\ntext that talks about the mechanisms in these lines..\nResponse: We acknowledge the importance of understanding the relationship between\nextreme heatwaves and low oxygen events, particularly in exploring the impact of\nheatwaves on oxygen solubility and stratification in waters.\nIn our revised manuscript, we have tried to address your concerns by investigating\nmodification of the oxygen solubility and stratification, thus elucidating their\nassociation with CLH events (Fig. R1). It can be found that the oxygen solubility in\nseawater during CLH events exhibits a significant decrease compared with it in the\naverage state. The regions with the most severe reductions in solubility are found in the\nequatorial Pacific and the mid-to-high latitude oceans. An enhanced stratification,\nrepresented by mixed layer shallowing, occurs in the widespread ocean during CLH\nevents. Most oceans experience a decrease in mixed layer depth, particularly\npronounced in the Southern Ocean and North Atlantic. Although some regions, such as\nthe equatorial Eastern Pacific, exhibit a slight increase in mixed layer depth, it remains\nnotably weak.\nFigure R1. The anomaly of oxygen solubility (a, unit: mmol/m3) and mixed layer depth\n(b, unit: m) during CLH events. Anomalies are obtained by subtracting the mean value\nfrom the value of these variables during the CLH event across each grid, for example,\nMLD = MLD - MLD . The seasonal variation cycle is removed before the\nanomaly CLH average\ncalculation, and results presented in this figure are derived from model simulations.\nIn conclusion, the decrease in oxygen solubility reduces the capacity for oxygen\nstorage in seawater. Meanwhile, ocean stratification intensifies during CLH events,\nwhich could hinder the ventilation of oxygen in the waters. The discussions mentioned\nabove can be found in the revised manuscript at Lines 319-326. However, due to time\nconstraints, we were unable to quantify the specific contributions of these factors. This\nwill be the focus of our future work.\nOnce again, we appreciate your thoughtful review and constructive feedback. If\nyou have any further suggestions or questions, please do not hesitate to reach out.\nSpecific Comments:\nEvery figure must be clear on the depth level that is being analsyed or compared.\nEvery figure in the responses to the authors never mentions the depth layer that is\nbeing shown. Depth is never mentioned in the responses text. There is only one place\nin the main text where the depth is revealed, and that is well into the methods section,\nand yet is should appear in every figure caption. Please introduce the surface depth\nrange and define the layer of interest very early on in the text, before line 64, and in\nthe abstract (page 1). Every figure caption should say surface (0-10 m). Much of the\noceanography community think about the full water column, whereas this paper is\nonly about the 0-10m surface layer of the ocean. Makes a huge difference in\ncommunication.\nResponse: Thanks for your suggestions. We have clarified in the abstract (Line 36) and\nat the beginning of the manuscript (Lines 90-91) that this work focuses on the surface\nlayer of the ocean (0-10m). Additionally, this information is stated in each figure (see\nFig. 1-4 in the revised manuscript).\nFigure R8. Once again this figure should say surface (0-10m) layer or similar. Is this\nthe trend from the beginning of the record (1960 to 2014) say? The panels of this\nfigure are showing that the correlation with the time trend in the data.\nResponse: In the revised manuscript, we have clarified that this figure pertains\nexclusively to the ocean’s surface layer (see Supplementary Fig. S20). It depicts the\ntemporal correlation coefficients between model simulations and observational data\nthroughout their overlapping periods, which span from 1982 to 2014 for temperature\nand from 2004 to 2014 for dissolved oxygen.\nI was and remain confused on what is presented in the figures and what is a result\nfrom a simulation, and what is a result from an observation? While it could easily be\nstated in the text, it helps enormously to know from the caption whether a particular\npanel is composed of observation an or simulations. See the questions about each\nfigure.\nResponse: We apologize for the lack of clarity regarding the source of the results\npresented in the figures. To rectify this issue, we have revised the captions to explicitly\ndistinguish between panels derived from observational data and those based on\nsimulation results. Please refer to our point-by-point responses for each figure below.\nFigure 1: This paper uses both models and observations. In this figure I cannot tell\nwhich panel is an observation? Which panel is a simulation. In panel d, I think this\nis simulation, but I cant see dots that are observations? I can see\nhistograms/ridgelines connected to dots on the map. Are histograms/ridgelines\ncomposed entirely of observations of temperature and oxygen in the regions whose\nboundaries are explained in the response to reviewers? Are they simulations?\nResponse: We acknowledge the importance of clearly differentiating between observed\ndata and simulated data in Fig, 1. We have updated the figure legend to explicitly label\neach panel, indicating whether it represents observations or simulations. Specifically:\nPanels a-c are derived from model simulations, and Panel d incorporates both\nobservational data and simulation results. The contour plot illustrating the global\ndistribution of LMF originates from simulations, while the associated histogram depicts\nLMF distribution across ocean grids in simulations. The remaining components of this\npanel, namely, the nine dots on the map and the interconnected ridgeline plots, are\nentirely based on observations of temperature and oxygen in the nine ocean basin\nregions.\nWe hope these changes address your concerns and improve the clarity of our\npresentation.\nFigure 2: This panel is clear that each panel is simulations. Results entirely from the\nensemble of 5 models.\nResponse: Thank you for your feedback on our manuscript. As you noted, this panel is\nderived from model simulations. We appreciate your attention to the detail.\nFigure 3: This is a complex figure. Are a-d observations or are the simulations? Or\nis there a mixture of simulations and observations in these panels? I can't tell from\nthe Figure caption. Same question for panels f-i as for a-d. Panel e: is this also\nsimulated fish biomass (eg from FISHMIP) or is it an observation? Panel J: not sure\nhow to interpret panel J, does it imply higher biomass is associated with longer CLH\ndays? While it shows big CLH are associated with high biomass, it doesnot show that\n6\nareas of high biomass are more affected by high CLH days than low biomass, and\nisnt this the more important result but shown in Fig S6? Compare with text in lines\n183:185Lines 192-193. Panel k: Why are there 6 regions (pie charts) when there are\nonly 4 biomass locations shown in panel e? Is Panel k observations or simulations?\nResponse: Panels a-d are derived from simulations, so are Panels f-i. The distribution\nof fish biomass in Panel e is sourced from Hatton et al. 2021, which are estimates based\non two fishery models, that is, BOATS (Carozza et al. 2017) and FEISTY (Petrik et al.\n2019). We have added an additional pelagic fish biomass data based on the global\nmarine fish model DBEM from a recently published article by Le Grix et al. (2023) in\nour first round of revision, and Fig. S6 is derived from it.\nIn Panel j, we categorized the ocean into several groups based on fish biomass\nlevels and then calculated the trend of CLH days within these groups. Results show that\nhigh fish biomass areas experience a more significant increase in the number of CLH\ndays compared to areas with lower biomass, suggesting that extreme events might have\na greater impact on regions with high biomass. Meanwhile, Fig. S6 demonstrates that\nfish biomass declines more rapidly with an increasing trend of CLH days. We believe\nthat Panel j and Fig. S6 illustrate the relationship between CLH events and fish biomass\nfrom different perspectives.\nSimilar to Panel j, in Panel k, we also divided the global ocean into six parts based\non biomass and calculated the proportion of each part where the trend of increasing\nCLH days is greater than the global average trend (0.96 days per year). This allows us\nto see whether the increase in extreme events on different biomass regions exceeds the\nglobal average level. The four locations selected in panel e are the iconic fishing\ngrounds in the ocean. Through the selection of these locations, we can observe the\ncharacteristics of changes in extreme events in these fishing grounds, which directly\nimpact the quantity of fish available for human consumption. In conclusion, Panel k\nclassifies the ocean into six sections based on marine biomass, while the four fishing\ngrounds selected in Panel e represent regions with abundant fish biomass,\ncorresponding to the final category in Panel k. All results in Panel k are derived from\nsimulations.\nReference\nCarozza D. A., Bianchi D. & Galbraith E. D. Formulation, general features and global\ncalibration of a bioenergetically-constrained fishery model. PLOS ONE 12,\ne0169763 (2017).\nHatton, I. A., Heneghan, R. F., Bar-On, Y. M. & Galbraith, E. D. The global ocean size\nspectrum from bacteria to whales. Sci. Adv. 7, 1–13 (2021).\nLe Grix, N., Cheung, W. L., Reygondeau, G., Zscheischler, J. & Frölicher, T. L. Extreme\n7\nand compound ocean events are key drivers of projected low pelagic fish biomass.\nGlob. Change Biol. gcb.16968 (2023).\nPetrik C. M., Stock C. A., Andersen K. H. & van Denderen P. D. Bottom-up drivers of\nglobal patterns of demersal, forage, and pelagic fishes. Prog. Oceanogr. 176,\n102124 (2019).\nFigure 4: Panel a: “anthropogenic component” is based on simulations (but through\nthe EEMD estimation. Panel b: is also based (and estimated from) on simulations. It\nstill needs stating. Units should be added to the colorbar. I like panel C. I am\nassuming this is from simulations (analysed 5 times from the 5 ensemble simulations)\nand the 9 modes climate are each diagnosed from each model simulation using the\nmodel indices? Is that right? I like the map in C and the bar chart. Is this map (and\nbar charts) estimated from the panel b natural variability or from the time series of\nboth panel a and panel b. What is the period for this map, I find it a bit odd to talk\nabout increased occurrence of CLH and connecting that directly to the dominant\nmodes. Does this include anthropogenic forcing or not? And for what time period is\nthe enhancement of occurrence estimated over, and for all of the ensemble runs?\nResponse: Thank you for your comments and suggestions. We have revised the figure\nto clarify that these results are based on simulations. Additionally, we have added units\nto the colorbar for better interpretation. True as you mentioned, in Panel c, the nine\nclimate modes are diagnosed from each model simulation using the model indices. The\nmap and bar chart are estimated from natural variability (shown in Panel b), which\nexclude the anthropogenic forcing component. The enhancement of occurrence is\nestimated over the period from 1985 to 2014 for all ensemble runs.\nThe authors have carefully revised the manuscript in response to the reviewers\ncomments and have added interesting and helpful analyses. In particular, I think the\nadditional supporting information analyses, such as comparing EEMD results to\nthose based on a shifting baseline, are very interesting. Overall, the methods\nunderlying the results are now much clearer to me. Despite these improvements, I\nstill have remarks and questions that should be addressed before publication.\nResponse: Thank you for your thoughtful and encouraging feedback on our manuscript.\nWe are pleased to hear that you found the additional analyses interesting and helpful,\nand that the methods underlying our results are now clearer. We appreciate your\ncontinued engagement with our work and the opportunity to further improve the\nmanuscript. Our point-by-point response to your remaining remarks and questions is\nprovided below.\nMajor comments:\nI am somewhat missing a more careful embedding of the present study in the\npublished research on low oxygen extremes and their interaction with marine\nheatwaves in the introduction section. I think this point is very important. Otherwise\nthe contribution of this paper, what it adds to existing knowledge on LOEs and CLH\nevents, can not be properly put into context. See more specific comments below.\nResponse: We appreciate your suggestion to highlight the need for a more thorough\nembedding of our study within the existing literature on low oxygen extremes and their\ninteraction with marine heatwaves. We agree that this context is crucial. According to\nyour suggestion, we have revised the Introduction section as follows:\n Expanded Literature Review\nWe have included a more comprehensive review of the current research. In the\nrevised manuscript, we highlight key findings and knowledge gap in the existing\nliterature.\n Specific Studies Referenced\nWe have added references to several seminal and recent studies in this research field.\nBased on your suggestions in specific comments, references such as Köhn et al. 2022\nand Mogen et al. 2022 have been included to embedding of the present study in the\npublished research on the low oxygen extremes and their interaction with marine\nheatwaves.\nAdditionally, to respond to your third major comment, we have included references\nto highlight the importance of studying extreme low-oxygen conditions in surface\nocean water (e.g. Brauko et al. 2020) as well as the sensitivity of marine biomass to the\neffects of heatwave and low oxygen (e.g. Tran and Johansen, 2023).\n9\n Contextualizing Our Contribution\nWe have clarified how our study builds on and differs from previous researches.\nSpecifically, we emphasize the innovation of our methodology and the new insights our\nresults provide into the compounding of heatwave and low-oxygen extreme events.\n Revised Objectives and Significance\nWe have revised objectives and significance sections to better articulate how our\nstudy addresses existing gaps in the literature and contributes to the broader\nunderstanding of these phenomena.\nThe Introduction section related to these changes has been revised as follows (you\ncan find in Lines 62-98 in the revised manuscript):\n“Characterized by the abrupt oxygen decline in the ocean, the low-oxygen extreme\nevent could cause species to reach their critical oxygen limits decades earlier than\nexpected from the deoxygenation (Heinze et al., 2023), eventually resulting in the loss\nof biodiversity and habitats (Köhn et al., 2022). The emergence of these low-oxygen\nextreme events is influenced by multiple factors, including anomalous water mass\ndistributions and ocean eddies etc. (Gruber et al., 2021; Atkins et al., 2022). Importantly,\ncase studies have suggested that marine heatwaves may act as catalysts for extreme\nlow-oxygen events: the abnormally low dissolved oxygen levels are found during a\nstrong marine heatwave known as “the Blob” (Mogen et al., 2022); In SW Atlantic, a\nmarine heatwave with sea surface temperatures over 30°C coincided with the regional\nfirst record of hypoxia (Brauko et al., 2020). Observations on the Baltic Sea reveal\nheatwaves can contribute to seawater oxygen deficiency, especially in the shallow areas\n(Safonova et al., 2024). These researches highlight the importance of exploring the\ncoupling between heatwaves and extreme low-oxygen events. Many marine ecosystems\nare highly vulnerable to multiple stressors and as a consequence may suffer more severe\ndamage if low-oxygen extreme and heatwave occur simultaneously or in close\nspatiotemporal proximity, as such joint occurrence has the potential to act\nsynergistically, leading to mutually reinforcing impacts on the marine organisms (\nBoyd\n& Brown, 2015; Laufkötter et al., 2017; Gruber et al., 2021). The combination of high\nwater temperature and low oxygen concentration has exerted influences on marine life\nsuch as fishes and scallops (Tomasetti et al., 2023; Tran & Johansen, 2023).\nNevertheless, so far we still lack a detailed knowledge about how heatwave and low-\noxygen extreme coupling across the ocean. A global assessment of the modification of\nthe compound events, as well as the drivers behind them, is urgently needed under\nongoing climate change.\nTo fill this gap, here we conduct a comprehensive analysis of the relationship\nbetween marine heatwaves and low-oxygen extreme events in the global ocean, using\nin-situ observations and a fully coupled Earth system model (Data descriptions and\n10\nvalidations of model could be found in “Methods”). The evolutions of the compound\nevents of the two extreme events are investigated in our study and a novel framework\nis established to identify the drivers contributing to their occurrences. We focus on\nextreme events within the uppermost ocean layer (0-10 meters) because conditions in\nthe surface ocean, such as temperature and oxygen levels, exert profound influences on\nmarine ecosystems, encompassing fish populations, as well as the marine biodiversity\n(Brauko et al. 2020; Mattiasen et al., 2020). Following the approach developed in ref.\n5,14,20, extreme events in temperature (oxygen concentration) are defined as a\nprolonged instance with anomalously warm (low-oxygen) water (See “Methods” for\nthe detailed definitions). We expect our findings to further deepen the understanding of\nthe synergy of extreme events in the oceans.”\nThank you once again for your valuable feedback. We hope that the revised\nIntroduction meets your expectations and provides a more comprehensive context for\nour study.\nReference\nAtkins, J., Andrews, O., & Frenger, I. (2022). Quantifying the Contribution of Ocean\nMesoscale Eddies to Low Oxygen Extreme Events. Geophysical Research Letters,\ne2022GL098672.\nBoyd, P. W., & Brown, C. J. (2015). Modes of interactions between environmental\ndrivers and marine biota. Frontiers in Marine Science, 2.\nBrauko, K. M., Cabral, A., Costa, N. V., Hayden, J., Dias, C. E. P., Leite, E. S., et al.\n(2020). Marine Heatwaves, Sewage and Eutrophication Combine to Trigger\nDeoxygenation and Biodiversity Loss: A SW Atlantic Case Study. Frontiers in\nMarine Science, 7, 590258.\nGruber, N., Boyd, P. W., Frölicher, T. L., & Vogt, M. (2021). Biogeochemical extremes\nand compound events in the ocean. Nature, 600(7889), 395–407.\nHeinze, C., Blenckner, T., Brown, P., Fröb, F., Morée, A., New, A. L., et al. (2023).\nReviews and syntheses: Abrupt ocean biogeochemical change under human-made\nclimatic forcing – warming, acidification, and deoxygenation.\nKöhn, E. E., Münnich, M., Vogt, M., Desmet, F., & Gruber, N. (2022). Strong Habitat\nCompression by Extreme Shoaling Events of Hypoxic Waters in the Eastern\nPacific. Journal of Geophysical Research: Oceans.\nLaufkötter, C., John, J. G., Stock, C. A., & Dunne, J. P. (2017). Temperature and oxygen\ndependence of the remineralization of organic matter. Global Biogeochemical\nCycles, 31(7), 1038–1050.\nMattiasen, E. G., Kashef, N. S., Stafford, D. M., Logan, C. A., Sogard, S. M., Bjorkstedt,\nE. P., & Hamilton, S. L. (2020). Effects of hypoxia on the behavior and physiology\n11\nof kelp forest fishes. Global Change Biology, 26(6), 3498–3511\nMogen, S. C., Lovenduski, N. S., Dallmann, A. R., Gregor, L., Sutton, A. J., Bograd, S.\nJ., et al. (2022). Ocean Biogeochemical Signatures of the North Pacific Blob.\nGeophysical Research Letters, 49(9), e2021GL096938.\nSafonova, K., Meier, H. E. M., & Gröger, M. (2024). Summer heatwaves on the Baltic\nSea seabed contribute to oxygen deficiency in shallow areas. Communications\nEarth & Environment, 5(1), 106.\nTomasetti, S. J., Hallinan, B. D., Tettelbach, S. T., Volkenborn, N., Doherty, O. W.,\nAllam, B., & Gobler, C. J. (2023). Warming and hypoxia reduce the performance\nand survival of northern bay scallops ( Argopecten irradians irradians ) amid a\nfishery collapse. Global Change Biology, 1–16.\nTran, L. L., & Johansen, J. L. (2023). Seasonal variability in resilience of a coral reef\nfish to marine heatwaves and hypoxia. Global Change Biology, 29(9), 2522–2535.\nI am a bit confused about the identification of correlations and definition of the\nextremes from the WOD data (e.g., Supplementary Figure S2) across the ocean\nbasins, or more general how the WOD data is handled. While you nicely explain how\nthe extremes are defined from the monthly data in the 9 basins, it is still illusive to\nme how the data is aggregated over the 9 basins in the first place. In the methods,\nyou refer to the Schmidtko study - however it is not explained yet what is done exactly.\nDo you just average all measurements in a month, say November 2008, over a whole\nbasin, to eventually get one monthly time series for each basin? If this is done, then\nit is questionable how this procedure influences correlations and compound event\noccurrence. For example it may be that in a November, only measurements from the\ngulf of Alaska enter the time series of the North Pacific, with much different mean\ntemperature and oxygen compared to the subtropical North Pacific. My thoughts\nhere may not apply, but a clarification is needed.\nResponse: We apologize that the current explanation in the methods section may lack\nsufficient detail, leading to potential confusion regarding how the WOD data is mapped.\nHere we’d like to further explain the data processing procedure. As you mentioned,\nwe referred to Schmidtko’s paper to complete a series of technical steps such as data\nfiltering and preprocessing (quality control, duplicate removal, etc.). As for how these\ndata were aggregated into 9 basins, we first divided the ocean into grids of 2.5 degrees\nby 2.5 degrees and mapped each profile data to the grid point closest to the recorded\nlatitude and longitude at the time of observation. Then, we followed the defined basin\nboundaries, performed a weighted average of all eligible data based on their grid\ncoordinates, and eventually obtained the monthly time series for each basin. As you\ncommented, this averaging is indeed inevitably affected by missing values. However,\n12\nlimited by the spatiotemporal coverage of WOD data, this is a compromise approach\nwe adopted after comprehensive consideration. Now we have added the descriptions of\nhow we mapped the WOD data and the limitations of this mapping to the methods\nsection (see Lines 354-363 in the revised manuscript). To ensure the reliability of our\nresults, we incorporated a machine learning-based dissolved oxygen dataset (GOBAI-\nO2) in the first round of revisions. The characteristics of dissolved oxygen from\nGOBAI-O2 show an overall consistency with our results (see model evaluation section\nin methods).\nLow oxygen extremes are often looked at in the ocean interior. I think the authors\nneed to establish why low oxygen extremes at the surface are relevant to make strong\nclaims such as in L71-L74. Can you give a reference for this claim? Likewise, the\nfinal sentence of the manuscript, 'Under their combined effects, marine heatwave\namplifies the peak signal of deoxygenation, highlighting urgent needs for adaptation\nresponses planning for future CLH risks.', is too strong if the importance of low\nsurface oxygen events can not be established.\nResponse: Thank you for your constructive comments. In response to your concern\nregarding the relevance of low oxygen extremes at the ocean surface (Lines 71-74), we\nagree that this aspect needs clearer justification. Typically, low oxygen extremes are\noften studied in the ocean interior. However, recent research has indicated that low\noxygen events in the near-surface ocean water are also critical, which could have\nsignificant ecological and biogeochemical implications and its relationship with marine\nheatwave should be particularly focused on. For example, a case study in the SW\nAtlantic reveals that a marine heatwave with sea surface temperatures over 30°C\ntriggered the region’s first recorded hypoxia event (Brauko et al., 2020). The synergistic\neffects of extreme temperatures and low oxygen affected microbenthic and\nphytoplanktonic biodiversity in this region. Mattiasen et al., 2020 demonstrates that the\nlow-oxygen extreme (measured in situ at 12 m depth) can have significant implications\non the rockfishes which live in shallow waters. Besides, Tran and Johansen, 2023\nshows the important impact of surface marine heatwaves and oxygen deficiency on fish\nresilience. These references support our claim about the relevance of studying these\nphenomena, particularly under the scenarios of exacerbated climate impacts. These\nreferences have been added to the Introduction in the revised manuscript, with details\nprovided in our responses to major comment 1.\nConcerning the final sentence of our manuscript, we acknowledge that the\nstatement might appear somewhat strong. In light of your feedback, we have revised\nthis sentence to better reflect the supported conclusions drawn from our data and to\nalign with the existing literature: “Our findings underscore the potential for heightened\n13\ncompound extreme risks to marine ecosystems, necessitating proactive planning in\nadaptation strategies for future climate-linked hazards.” (Lines 340-342 in the revised\nmanuscript)\nAdditionally, we plan to investigate the impacts of marine heatwaves on extreme\nlow-oxygen events in the ocean interior in our future research. It is interesting to see\nthat enhanced stratifications, indicated by mixed layer shallowing, occurs in the ocean\nduring the events (Fig. R2), which implies that oxygen in the ocean interior might be\naffected as a consequence of suppressed oxygen supply under heatwaves. However,\ninvestigating these effects is difficult due to the diminishing availability of\nobservational data at greater depths. Our future research will concentrate on addressing\nthese issues.\nFigure R2. The anomaly of mixed layer depth (unit: m) during CLH events. Anomalies\nare obtained by subtracting the mean value from the value of the variable during the\nCLH event (i.e., MLD = MLD - MLD ) across each grid. The seasonal\nanomaly CLH average\nvariation cycle is removed before the calculation, and results presented in this figure\nare derived from model simulations.\nReference\nBrauko, K. M., Cabral, A., Costa, N. V., Hayden, J., Dias, C. E. P., Leite, E. S., et al.\n(2020). Marine Heatwaves, Sewage and Eutrophication Combine to Trigger\nDeoxygenation and Biodiversity Loss: A SW Atlantic Case Study. Frontiers in\nMarine Science, 7, 590258.\nMattiasen, E. G., Kashef, N. S., Stafford, D. M., Logan, C. A., Sogard, S. M., Bjorkstedt,\nE. P., & Hamilton, S. L. (2020). Effects of hypoxia on the behavior and physiology\nof kelp forest fishes. Global Change Biology, 26(6), 3498–3511\nTran, L. L., & Johansen, J. L. (2023). Seasonal variability in resilience of a coral reef\nfish to marine heatwaves and hypoxia. Global Change Biology, 29(9), 2522–2535.\nSpecific comments:\nsee, only analyze one climate model simulation, this should be \"and a state-of-the-\nart climate model simulation\"\nResponse: Thank you for your careful review. We have revised the text and this phrase\nhas been changed to \"and a state-of-the-art climate model simulation.\"\nL39: \"Our work shows a widespread emergence of low-oxygen extreme events in\nresponse to increasing marine heatwaves\" -> as you write in the response to\nreviewers, and have already corrected earlier in the abstract, isn't it rather mainly an\nincrease in low-oxygen extreme events due to increasing mean temperatures than\ndue to an increase in marine heatwaves?\nResponse: We appreciate your thoughtful comments on the statement here. In response\nto your feedback, we have revised the abstract to more accurately reflect the findings\nof our study. The revised statement (Lines 40-42) as follows:\n“The rise in the compound events is primarily attributable to long-term warming\nprimarily induced by anthropogenic forcing, in tandem with natural internal variability\nmodulating their spatial distribution.”\nL58: I would propose adding Hobday et al., 2016 and Frölicher et al 2018 to the\nreferences on marine heatwaves\nResponse: Following your suggestions, we have added these references in the revised\nmanuscript.\nL58-59: There is some work on low oxygen events (as you also write before) and also\nto a lesser extent on CLH events, for example Köhn et al on oxygen extremes, this\nshould be added here as an example of studying low oxygen extremes, and also\nGruber et al. 2021, as well as Mogen et al 2022\n(https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2021GL096938) who\nspecifically looked at the heatwave impact on oxygen. Overall, I think the results of\nthe current study needs to be more carefully embedded in existing literature. I would\npropose going through section 5.3 of this recent discussion paper\nhttps://bg.copernicus.org/preprints/bg-2023-182/ for further references for studies\nabout low oxygen extremes. The knowledge gap needs to be stated more precisely.\nResponse: Thank you for your detailed comments and for suggesting several key\nreferences that add substantial depth to our discussion of low oxygen and CLH events.\nWe appreciate your guidance in enhancing the contextual foundation of our study within\nthe existing literature. To address your concern, we have revised this section and added\nthese references to provide a more comprehensive background and support for our\n15\nstudy’s context, with details provided in our responses to major comment 1. We hope\nthese revisions adequately could respond to your suggestions and improve the\nmanuscript.\nFigure 1: The corresponding correlation coefficients are presented in the top right\ncorner with symbols * and ** representing the statistically significant at P<0.01 and\nP<0.1 level, respectively.\n-> * and ** were switched in this sentence, ** should refer to P<0.01 I think.\nResponse: Thank you for pointing out the error here. We have corrected it in the revised\nmanuscript.\nL79-L82: The additional analysis on the correlation coefficient and LMF if the data\nis detrended (Text S1 and Figure S21) provides very important information. Later on,\nyou discuss that the LMFs > 1 are weaker but still present when detrending the data\n(L135-139). I would propose to add also somewhere that the correlations of T and O2\nanomalies (shown in Fig 1d) are less negative but still present after detrending the\ndata.\nResponse: In response to your suggestion, we have added this statement in Lines 103-\n105 in the revised manuscript:\n“The correlation between temperature and oxygen concentration anomalies becomes\nless negative yet remains evident after data detrending, as detailed in Supplementary\nText S1.”\nL83-86: Thank you for looking into the other metrics in response to my comments. I\nam not yet fully convinced by the reply in the R2R about why the event frequency is\nthe best choice here. Given that you show a composite of the number of extreme days\nfor single and overlapping events in the LMF map in Fig 1d, the most useful choice\nfor Fig 1a,b would, in my opinion, be to also show number of event days here. The\nfrequency of individual events is not only influenced by how frequently the variables\nare above the thresholds (i.e., the number of days), but also inversely controlled by\nhow long individual events take on average, and thus somewhat harder to interpret.\nResponse: Thank you for your thoughtful comments on the choice of metrics here. We\nconcur with your suggestion to focus on the number of event days in Fig. 1a and 1b,\nrather than on the frequency of events. Accordingly, we have updated these figures in\nthe revised manuscript to illustrate the number of extreme event days, as shown in Fig.\nR3 here. Figures that previously displayed the frequency of events have been relocated\nto Fig. S3 in the supplementary material.\n1999 (a) and 2000-2014 (b). For sake of comparisons, the number of days for each\nextreme event type is classified into four levels, that is, Rare (R), Occasional (O),\nFrequent (F), Common (C), respectively.\nL93-95: Fig 1a,b is clearly interesting, but it does not directly provide proof about the\nclaim you make here. You show that areas with increasing MHW number are often\nalso areas with increasing LOE number. However, this does, strictly speaking, not\nimply that MHWs and LOEs co-occur more often, just that both events individually\nbecome more frequent in a region. You discuss co-occurrence starting L107ff and\nFig. 2 - I would leave it out here where it is not yet shown.\nResponse: Thank you for your insightful comments concerning the interpretation of\nthe data here. We fully agree with your suggestions and have accordingly removed these\nsentences from this section.\nL121-122: How is the global LMF estimated from the WOD data - is this the area-\nweighted average of the colored dots in Fig 1? and what is the global mean LMF\nfrom the CESM model? this information could be provided in L129\nResponse: Yes, the global LMF is estimated from area-weighted average of the data.\nThe global mean LMF from model simulation, according to our calculation, is\napproximate 5.1. We have updated this information in the revised manuscript (see Lines\n148-151):\n“Model simulations are in good agreement with the observations, showing that over 90%\nof the ocean experience more CLH events (i.e., LMF >1, contour lines in Fig. 1d), with\na global average LMF value of approximately 5.1.”\nFigure 1: It would be nice if the maps could be all centered the same way. This would\nfor example help identifying that the regions with high number of CLH days (Fig 2a)\nare also those with high LMF (Fig 1d).\nResponse: We appreciate your attention to these details and your help in improving the\npresentation and clarity of our work. We have revised Figure 1 to ensure all maps are\nuniformly centered.\nand which are from CESM.\nResponse: Thank you for your suggestion. Results in Fig. 2 are all derived from model\nsimulation, and we have clarified the data source in the caption:\n“Figure 2. Annual event days for compound low-oxygen extreme and heatwave (CLH),\nmarine heatwave only (MHW-only) and low-oxygen extreme only (LOE-only) in ocean\nsurface layer based on the model simulations.”\nL144: what is considerable increase? this should be specific when also given a specifc\nocean area (> 70%) for it.\nResponse: In response to your feedback, we have revised the text to specify the increase\nmore clearly. Here is the revised text (see Lines 164-167 in the revised manuscript):\n“Comparisons between two 15-year periods (1985-1999 and 2000-2014) in the climate\nmodel also show a considerable rise in the annual CLH days across more than 70% of\nthe global oceans, with an average increase exceeding 20 days, despite large spatial\nvariations exist.”\nL161-164: I see that the CLH fraction of LOEs is increasing, however, the proportion\nof marine heatwaves that are also LOE does not seem to increase at the same time\n(47% at the beginning and at the end). This somewhat weakens the statement made\nhere in my opinion.\nResponse: Thank you for your careful review and insightful comments. We calculated\nthe CLH fraction of MHWs during the three periods, which are 46.8%, 47.5% and 47.5%\nrespectively. As you noted, this increase is relatively modest and not as substantial as\nthe rise in the CLH fraction of LOEs.\nFrom our perspective, concerning low-oxygen extreme events, the degree to which\nthey are affected by heatwaves is steadily increasing. This is reflected in the Venn\ndiagram in Fig. 2b, where the blue circle continuously approaches and overlaps with\nthe red circle over time (or one could say the red circle gradually “engulfs” the blue\ncircle). Meanwhile, the red circle itself is also enlarging rapidly, resulting in a relatively\nmodest increase in CLH fraction of MHWs, which is the phenomenon you mentioned.\nIn conclusion, we think that the overall coupling degree of these two extreme events\nhas increased, albeit asymmetrically in heatwaves and extreme hypoxic events.\nHere we present an example to illustrate this situation. Suppose there is a period\nwhen the number of days with heatwaves and low-oxygen extreme events is 44 and 42\ndays, respectively, and their concurrent days (days of CLH events) total 20 days. This\ncan be represented by the Venn diagram in Fig. R4a: the proportion of CLH days to\n18\nheatwave days is 45.5% (20/44), and to low-oxygen extreme days is 47.6% (20/42).\nAfter some time, the number of heatwave days increased to 86 days, and low-oxygen\nextreme days to 51 days, with their concurrent occurrence, the CLH days, increasing to\n41 days, as shown in Fig. R4b. At this point, the proportion of CLH days to heatwave\ndays was 47.7% (41/86), and to low-oxygen extreme days was 80.4% (41/51). Although\nthe CLH fraction of MHWs does not increase substantially, the coupling degree of low-\noxygen extremes with heatwaves in Fig. R4b is obviously higher than in the previous\nFig. R4a.\nFurthermore, we’d like to emphasize that heatwaves can promote occurrence of\nlow-oxygen extreme events, but the impact of low-oxygen extreme events on the\noccurrence of heatwaves is minimal. Hence, from a physical significance perspective,\nwe think we should put our focus on the proportion of CLH events to low-oxygen\nextreme events, as it reveals the extent to which degree the heatwave influence low-\noxygen extreme events, rather than the proportion of CLH events to heatwaves.\nFigure R4. An example to illustrate coupling degree of low-oxygen extreme with\nheatwave. The blue (red) circle represents low-oxygen extreme events (marine\nheatwave).\nWe acknowledge that the expressions here are indeed somewhat coarse. The\ncorresponding sentences have been revised in the manuscript to describe our findings\nmore precisely (Lines 182-183 and 190-192):\n“the low-oxygen extreme events are becoming increasingly associated with heatwaves\nover the past few decades: …”\n“the joint occurrences (i.e., CLH events) have played a majority role in timeline of low-\noxygen extreme events, reflecting a substantial increase in its synchronization with\nheatwaves.”\nWe hope this response could adequately address your concerns and thank you once\nagain for helping us improve the precision and clarity of our analysis.\nL180: the normalization in Fig 3b,d,g,i - are the two periods normalized\nindependently from each other? this is relevant to interpret the percentage changes\nin the normalized anomalies.\nResponse: Thank you for your query regarding the normalization process for the data\npresented in Fig. 3. To clarify, the two periods were not normalized independently from\neach other in our analysis. Instead, we applied a consistent normalization approach for\nall values (whole period) to ensure the comparability of anomalies. Subsequently, the\nnormalized values were colored according to the period (before 2000 or after 2000)\nthey belong to.\nL185: \"reaching a peak of 2.53 days per year. This rise is approximately 2.6 times\ngreater\" I think it would be more informative to compare mean CLH trends over high\nfish biomass regions with global mean treands, rather than comparing a peak to a\nthe overall mean.\nResponse: We fully agree with your suggestion regarding the comparison of mean CLH\ntrends over high fish biomass regions with global mean trends rather than focusing on\na peak value compared to the overall mean. According to our calculation, the mean\nCLH trend over high fish biomass regions (>50 gm-2) is about 1.47 days per year. We\nhave added this result into the reivsed manuscript (see Lines 208-210):\n“In regions where fish biomass exceeds 50 gm-2, there is an upward trend of CLH days\nabout 1.47 days per year. This rise is approximately 1.5 times greater than the global\naverage increase.”\nL223f: The decomposition of the changes in CLH events shown in Fig 4a,b is very\ninteresting. However, while it is clear to me how Fig 4b is calculated (remove the\nanthropogenic forcing signal from the time series using the EEMD approach), I do\nnet yet understand how Fig 4a is calculated exactly. Some more explicit information\nwould be good here. Maybe as a sentence appended to the methods section paragraph\non EEMD?\nResponse: Thank you for your comments on the decomposition analysis presented in\nFig. 4. In response to your suggestion, we agree that additional information would\nbenefit the understanding of the methods employed.\nFig. 4a illustrates the influence of anthropogenic forcing on changes in annual CLH\ndays. Using the EEMD method, the original time series is decomposed into two\ncomponents: the anthropogenic forcing signal and the internal natural variability signal.\nBy conducting extreme event detection on the former, we can identify extreme events\nrelated to anthropogenic forcing. Figure 4a is derived by calculating the difference\nbetween the results from the two periods.\nWe apologize for the ambiguous description of the method for calculating extreme\nevents related to anthropogenic forcing. Relevant descriptions have been added to the\nmethods section, after the paragraph discussing the EEMD approach (see Lines 545-\n550):\n“Using the EEMD method, the original time series is decomposed into two components:\nthe anthropogenic forcing signal and the internal natural variability signal. By\nperforming extreme event detection on the anthropogenic forcing signal, we can\nidentify extreme events related to anthropogenic factors. Similarly, extreme events\nrelated to natural internal variability are identified by conducting extreme event\ndetection on the latter.”\nL247: 2000 or 1980? since referring to Fig 4b, it should be 2000 I think. Furthermore\nplease add the time period analyzed for Fig 4c to the caption.\nResponse: We appreciate your identification of our error here; indeed, as you noted, it\nshould be 2000. The time period analyzed for Figure 4c is 1985-2014. We have added\nthis information to the caption.\nL261f: a smaller thing: what is meant with \"not significant compared with other\nclimate modes\"? Usually statistical significance would be assessed for each mode\ncontribution individually I think\nResponse: We apologize for any misunderstandings caused by our unclear descriptions.\nOur intended meaning here is that the influence of this climate mode is relatively minor\ncompared to other modes, and therefore, it is not depicted in the figure. We have revised\nthis sentences which you can see in Lines 283-284:\n“However, the influence of NAO is relatively minor compared to other modes;\nhence\nit is not drawn in the figure.”\nL267: I would maybe also add that the former is a main contributor at the global\nscale (compared to the internal variability that is more important for the spatial\npatterns of changes)\nResponse: In response to your suggestion, we have revised the text to explicitly state\nimportant contribution of anthropogenic forcing at the global scale. The revised\nsentence is as follows:\n“Our results reveal that both anthropogenic forcing and natural internal variability are\nresponsible for the modifications of CLH events: the former is the main contributor to\nthe lengthening of annual CLH days at the global scale, while the latter plays an\nimportant role in spatial distributions of this change.”\nL275f: I agree that air-sea oxygen is most likely very important - however I think the\nexplanation here may not be very clear yet. I assume the variations in surface oxygen\ncontent from air-sea oxygen exchange or not so much driven by variations in\natmospheric oxygen but by changes in oxygen solubility, mainly due to changes in\nSST, causing a change in the capacity of the surface water to hold oxygen.\nResponse: Thank you for your valuable comments. We fully agree with your viewpoint\nthat changes in oxygen solubility due to sea temperature variations plays a crucial role\nin modifying surface ocean oxygen inventories. We apologize for overlooking this\naspect in our discussions. In response to your comments, we have revised the related\nparagraph as follows (Lines 304-309 in the revised manuscript):\n“For the oxygen concentration, changes in oxygen solubility due to sea surface\ntemperature variations plays a crucial role in modifying the seawater’s capacity to store\noxygen. Concurrently, the air-sea oxygen flux, along with the horizontal and vertical\ntransport of dissolved oxygen as well as local processes of oxygen production and\nconsumption, collectively shape surface dissolved oxygen concentrations.”\nFurthermore, we calculated the anomaly of oxygen solubility during the events (Fig.\nR5). It can be found that the oxygen solubility in seawater during CLH events exhibits\na significant decrease compared with the average state. The regions with the most\nsevere reductions in solubility are found in the equatorial Pacific and the mid-to-high\nlatitude oceans. The related discussions are added in Lines 319-323 in the manuscript.\nFigure R5. The anomaly of oxygen solubility (unit: mmol/m3) during CLH events.\nSimilar to Fig. R2, anomalies are obtained by subtracting the mean value from the value\nof the variable during the CLH event across each grid.\nL287: the suppression of oxygen supply through stratification applies to subsurface\nwaters, but I don't think it is important for surface waters that are studied here. Please\nrephrase.\nResponse: Thank you for your insightful comments regarding our discussion on the\nimpact of ocean stratification. We totally agree with your opinion that stratification\nprimarily influences oxygen levels in subsurface waters. However, other reviewers\nwant us to see how ocean stratification changes during the event. Therefore, in the\nrevised manuscript, we rephrased this section and transformed it as an additional\ndiscussion on how heatwaves might affect dissolved oxygen in the subsurface ocean\nthrough stratification (Lines 323-326 in the revised manuscript).\n“Interestingly, significant enhanced stratification, indicated by mixed layer shallowing,\noccurs in the ocean during CLH events (Supplementary Fig. S11b, referred to as Fig.\nR2 in our response). This intensified ocean stratification can hinder the ventilation of\noxygen, suppress the ocean’s oxygen supply, which might eventually affect oxygen in\nthe ocean interior.”\nL292-295: Do you show that the changes during the extremes are larger than those\nexpected from decades of anthropogenic changes? This statement may be a bit strong\nand not hold generally - please check.\nResponse: Thank you for your valuable comments. We have carefully reviewed the\nrelevant sentences and acknowledge that, as you pointed out, the original statement here\nis overly strong. Therefore, we have revised it to better align with our paper:\n“the associated extreme events can deteriorate the living environment within a\nrelatively short period of time, resulting in devastating impacts on marine ecosystems\nand biota.”\nL295f: You show some trend that LOEs tend to be more often also MHWs over time,\nbut it does not seem like MHWs are more often also LOEs (my comment above),\n'intensified synergy' is maybe thus a bit too coarse?\nResponse: As mentioned in our reply to your comments on L161-164, we acknowledge\nthat there is only a modest increase in the CLH fraction of heatwave days from 46.8%\nto 47.5%. The expressions here are indeed somewhat coarse. To describe the relevant\nresults more precisely, we revised the sentence to: “The low-oxygen extreme events are\nbecoming increasingly associated with heatwaves.”\nSupporting information text S3: This addition is very useful. I just want to make one\ncomment as just my perspective on it. A shifting baseline, in my understanding, does\nnot necessarily need to be based on a shifting window for extreme event definition\n(that is discussed in this section). It can also rather be a shifting-mean baseline -\nalready realized by simply detrending the data based on some linear or polynomial\nspline, or subtracting the ensemble mean to identify forced changes (e.g., Burger et\nal, 2022), some low-pass like filtering, or the EEMD approach you use. After all, I\nwould consider the EEMD a (notably very interesting) option to identify the forced\nchanges for applying a shifting-mean baseline.\nResponse: Thank you for your constructive feedback on our discussion of shifting\nbaselines in the context of defining extreme events. We appreciate your detailed\nunderstanding and the alternative methods you provided.\nWe recognize the value in discussing alternative approaches to provide a\ncomprehensive view. To this end, we have expanded our discussion in Supporting\nInformation Text S3 to include a brief overview of the methods you mentioned. This\nwill not only enrich our manuscript but also help readers appreciate the range of\ntechniques available for handling similar analyses:\n“There are several methods for separating relatively short-lived and transient changes\nfrom long-term signals, such as detrending the data using linear or polynomial splines,\nsubtracting the ensemble means, and applying low-pass filtering. Here, we compare the\ncharacteristics of marine heatwaves derived from the shifting baseline approach, which\nuses a shifting window for extreme event definition, with those derived from the\nEEMD-based method.”\nThank you again for your insightful suggestions that will undoubtedly improve the\nmanuscript.\nMinor comments:\nL79 and at other places: The authors write \"extreme event\", although it should be\n\"extreme events\" in the plural form, if I am not mistaken.\nResponse: Thank you for pointing out the grammatical inconsistency related to the use\nof “extreme event” versus “extreme events”. We have corrected this error at line 79 and\nall other places it occurred. The revised manuscript now consistently uses “extreme\nevents”.\nL82: \"Outputs from climate model\" -> from A climate model\nResponse: Done.\nL167: Typo \"although slightly difference exists in the magnitudes\"\nResponse: Done.\nL189: \"(approximately +0.96 days per year)\" This information is already given\nbefore.\nResponse: Thank you for your comments. This repetition has been removed in the\nrevised manuscript.\nReviewer #1 (Remarks to the Author):\nI have re-read the manuscript and assessed the changes (as discussed in the response to\nreviewers). The authors have responded to my second review comments. I found the manuscript in\nvery good shape and I am satisfied with this revised version. The captions to each of the figures are\nclear. The distinction between trends and attribution to climate change is now accurate and clear.\nThe manuscript reads very well, no further changes from me.\nReviewer #2 (Remarks to the Author):\nThe authors have very carefully responded to all my comments. Thank you! I have one final science\nrelated question, probably beyond the scope of the article: Why is the rise in CHL events larger in\nhigh-biomass regions? It would be very interesting to provide some interpretation explaining this\nfinding.\nI support publication of the manuscript as is. However, reading it once more, I have a couple of\nsmall text-related comments:\n* L209: \"This rise is approximately 1.5 times greater than the global average increase.\" I find 1.5\ntimes greater to be somewhat ambiguous. I would propose to rephrase \"The rise is approximately\n50% greater than the global average increase.\"\n* L59: \"focusing on extreme warming event\" -> plural events. At some places, the choice of plural\nversus singular forms should be checked before final publication, but that may be also something\nfor the proofs.\n* L104: 'yet remain evident'\n* L111ff: \"the pattern of low-oxygen extreme events broadly resembles that of heatwave events,\nwith a pattern correlation coefficient over 0.8 (p < 0.01).\" Maybe it would be good to be more\nspecific here, writing something like 'the pattern of low-oxygen extreme event frequency and\nduration broadly resembles that of marine heatwave events, with pattern correlation coefficients\nover 0.8' (I presume it is two correlation coefficients)\n* L132: 'at each grid' do you mean something like 'in each grid cell'?\n* Fig S4: maybe you want to add to the caption of Fig S4 that panel b is the same as Fig 2a in the\nmain text, may facilitate building bridges to the main text.\n* L165f: \"show a considerable rise in the annual CLH days across more than 70% of the world’s\noceans, with an average increase exceeding 20 days\" do the 70% refer to the fraction of the ocean\nthat experiences any increase in CHL events? proposal, sth like: \"show a considerable rise in\nannual CLH days, with increases in more than 70% of the world's ocean and a global average\nincrease exceeding 20 days.\"\n* L171: something to double check: how can the global increase over 15 years be larger than 20\ndays (L165f), while the linear trend is about 1 per year globally, which would propose 15 days\nincrease in 15 years.\n* L188 'the strengthening'\n* L190-192: not sure whether I understand 'a majority role in timeline of low-oxygen extreme events',\nmaybe rephrase\n* L217-220: I am not sure I understand the sentence starting with 'This asymmetry..'. I would\npropose to rewrite the sentence. The finding of increased fish biomass loss where also CLH event\ndays increase more is very interesting. It is clear that no causal relationship \"CLH event increase ->\ndecline in fish biomass\" is necessarily implied, but it may open interesting scientific questions for\nfuture research.\n* L283f: to clarify \"However, the influence of NAO is relatively minor compared to other modes;\nhence it is not drawn in the figures.\" NAO is shown on the map in panel c, it's just not represented by\na separate bar plot subpanel, right?\n* L284ff: I struggle to match this sentence with Fig. S8a, there does not seem to be an influence\nclose to Greenland. Please double-check\n* L297: \"though we acknowledge that some signals (such as aerosols and volcanic activities) are\nstill mixed in\" lacks some clarity, are you expressing they influence both anthropogenic and internal\nvariability components? maybe rephrase\n* L334: Should be \"a marked increase in\", I think.\n* L433-434: \"Nevertheless, the simulated distribution is in good agreement with it in observations\nfrom WOA 2018\" WOA or WOD? in the beginning of the section you refer to WOD. To me it is not\nclear whether you compare the model against WOD or WOA (which is of course mainly building on\nWOD) in this section, or compare to both as indicated by the section title?\nResponse to Reviewer #1’s Comments\nI have re-read the manuscript and assessed the changes (as discussed in the response\nto reviewers). The authors have responded to my second review comments. I found\nthe manuscript in very good shape and I am satisfied with this revised version. The\ncaptions to each of the figures are clear. The distinction between trends and\nattribution to climate change is now accurate and clear. The manuscript reads very\nwell, no further changes from me.\nResponse: Thank you for your thorough re-evaluation of our manuscript and for your\npositive feedback on the revisions we have made. We are pleased to hear that you found\nthe manuscript in very good shape and overall readability have met your expectations.\nWe greatly appreciate your constructive comments throughout the review process,\nwhich have been instrumental in improving the quality and clarity of our work.\nResponse to Reviewer #2’s Comments\nThe authors have very carefully responded to all my comments. Thank you! I have\none final science related question, probably beyond the scope of the article: Why is\nthe rise in CHL events larger in high-biomass regions? It would be very interesting\nto provide some interpretation explaining this finding.\nResponse: Thank you for your positive feedback and for raising an interesting scientific\nquestion. We think that the larger rise in CLH events in these regions can be attributed\nto several factors. High-biomass areas are primarily located in near-coastal regions.\nThese areas typically have more abundant nutrient availability, which can support more\nfrequent and intense phytoplankton blooms, including those associated with CLH\nevents. Moreover, Guo et al. (2022) revealed that due to boundary currents, coastal\nprocesses, and ocean eddy fluxes, the heatwaves experienced in high-biomass regions\nmay be more severe than expected, which also provides a possible explanation for your\nquestions.\nWe agree that providing a detailed interpretation of this finding would be very\ninteresting and could be a valuable direction for future research. Although this is beyond\nthe scope of our current manuscript, we appreciate the suggestion and will consider it\nfor further investigation.\nThank you again for your insightful comments and contributions to our work.\nReference\nGuo, X., Gao, Y., Zhang, S., Wu, L., Chang, P., Cai, W., et al. (2022). Threat by marine\nheatwaves to adaptive large marine ecosystems in an eddy-resolving model. Nature\nClimate Change, 12(2), 179–186.\nI support publication of the manuscript as is. However, reading it once more, I have\na couple of small text-related comments:\nResponse: Thank you for your support for the publication of our manuscript. We\nappreciate your attention to detail and are grateful for your feedback. We have carefully\nrevised our manuscript and addressed the comments you have provided below to\nimprove the clarity and quality of the manuscript.\n* L209: \"This rise is approximately 1.5 times greater than the global average\nincrease.\" I find 1.5 times greater to be somewhat ambiguous. I would propose to\nrephrase \"The rise is approximately 50% greater than the global average increase.\"\nResponse: The sentence has been revised to: \"The rise is approximately 50% greater\nthan the global average increase.\"\n* L59: \"focusing on extreme warming event\" -> plural events. At some places, the\nchoice of plural versus singular forms should be checked before final publication,\nbut that may be also something for the proofs.\nResponse: Done.\n* L104: 'yet remain evident'\nResponse: Done.\n* L111ff: \"the pattern of low-oxygen extreme events broadly resembles that of\nheatwave events, with a pattern correlation coefficient over 0.8 (p < 0.01).\" Maybe it\nwould be good to be more specific here, writing something like 'the pattern of low-\noxygen extreme event frequency and duration broadly resembles that of marine\nheatwave events, with pattern correlation coefficients over 0.8' (I presume it is two\ncorrelation coefficients)\nResponse: We have revised this sentence following your suggestion.\n* L132: 'at each grid' do you mean something like 'in each grid cell'?\nResponse: Yes, it should be ‘in each grid cell’, we have corrected it in the revised\nmanuscript.\n* Fig S4: maybe you want to add to the caption of Fig S4 that panel b is the same as\nFig 2a in the main text, may facilitate building bridges to the main text.\nResponse: Thank you for your valuable suggestion. We have added this information to\nthe caption of Fig. S4.\n* L165f: \"show a considerable rise in the annual CLH days across more than 70%\nof the world’s oceans, with an average increase exceeding 20 days\" do the 70%\nrefer to the fraction of the ocean that experiences any increase in CHL events?\nproposal, sth like: \"show a considerable rise in annual CLH days, with increases in\nmore than 70% of the world's ocean and a global average increase exceeding 20 days.\"\nResponse: Yes, the 70% refer to the fraction of the ocean that experiences any increase\nin CLH events. Following your suggestion, this sentence has been revised to “show a\nconsiderable rise in annual CLH days, with increases in more than 70% of the world's\nocean and a global average increase exceeding 20 days.”\n* L171: something to double check: how can the global increase over 15 years be\nlarger than 20 days (L165f), while the linear trend is about 1 per year globally, which\nwould propose 15 days increase in 15 years.\nResponse: We have carefully checked the statement here. In line 165, the change is\ncalculated based on the difference between the average of two 15-yr period 1985-1999\nand 2000-2014 (that is, avg minus avg ), which is therefore exceeding\n_2000-2014 _1985-1999\n20 days. This does not contradict the result (linear trend of 0.96 day per year) described\nin L171.\n* L188 'the strengthening'\nResponse: Done.\n* L190-192: not sure whether I understand 'a majority role in timeline of low-oxygen\nextreme events', maybe rephrase\nResponse: We apologize for any confusion caused by the original phrasing. To improve\nclarity, we have rephrased the sentence as follows: “The joint occurrences (i.e., CLH\nevents) has become the dominant concern in the evolution of low-oxygen extreme\nevents”\n* L217-220: I am not sure I understand the sentence starting with 'This asymmetry..'.\nI would propose to rewrite the sentence. The finding of increased fish biomass loss\nwhere also CLH event days increase more is very interesting. It is clear that no causal\nrelationship \"CLH event increase -> decline in fish biomass\" is necessarily implied,\nbut it may open interesting scientific questions for future research.\nResponse: The term “asymmetry” here refers to the greater increase of extreme events\nin high-biomass areas compared to low-biomass areas, illustrating an asymmetrical\nincrease pattern between these regions. To enhance clarity, we have rewritten the\nsentence as follows: “This greater increase of CLH events in high-biomass regions\nimplies more severe damage to marine life and ecosystem”. Our discussions about why\nthe rise in CLH events is larger in high-biomass regions can be found in our response\nto Comment #1.\n* L283f: to clarify \"However, the influence of NAO is relatively minor compared to\nother modes; hence it is not drawn in the figures.\" NAO is shown on the map in panel\nc, it's just not represented by a separate bar plot subpanel, right?\nResponse: Yes, you are totally right.\n* L284ff: I struggle to match this sentence with Fig. S8a, there does not seem to be\nan influence close to Greenland. Please double-check\nResponse: We have carefully checked this sentence. The NAO indeed has some impact\non the oceans east of Greenland, as shown in Fig. R1 (despite the affected area is\nrelatively limited). Nevertheless, to avoid any confusion, we have removed this\nsentence.\nFigure R1. Percentage change in annual CLH event days in NAO positive phase. This\nfigure is same as Fig. S8a in the Supplementary material.\n* L297: \"though we acknowledge that some signals (such as aerosols and volcanic\nactivities) are still mixed in\" lacks some clarity, are you expressing they influence\nboth anthropogenic and internal variability components? maybe rephrase\nResponse: We apologize for the lack of clarity in this sentence. Here we intend to\nindicate that our method exclusively separates extreme events into anthropogenic and\ninternal variability components; it does not identify specific signals such as aerosols\nand volcanic activities. We have rephrased this sentence as follows: “our method does\nnot identify specific signals such as aerosols and volcanic activities”.\n* L334: Should be \"a marked increase in\", I think.\nResponse: Thank you for your comment. We have revised this sentence.\n* L433-434: \"Nevertheless, the simulated distribution is in good agreement with it in\nobservations from WOA 2018\" WOA or WOD? in the beginning of the section you\nrefer to WOD. To me it is not clear whether you compare the model against WOD or\nWOA (which is of course mainly building on WOD) in this section, or compare to\nboth as indicated by the section title?\nResponse: As described in Line 343 and 362, both WOD and WOA are used in model\nvalidation. We compare the model results against to both of them as indicated by the\nsection title."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I have re-read the manuscript and assessed the changes (as discussed in the response to reviewers). The authors have responded to my second review comments. I found the manuscript in very good shape and I am satisfied with this revised version. The captions to each of the figures are clear. The distinction between trends and attribution to climate change is now accurate and clear. The manuscript reads very well, no further changes from me.", "reviewer_2": "The authors have very carefully responded to all my comments. Thank you! I have one final science related question, probably beyond the scope of the article: Why is the rise in CHL events larger in high-biomass regions? It would be very interesting to provide some interpretation explaining this finding.\nI support publication of the manuscript as is. However, reading it once more, I have a couple of small text-related comments:\n* L209: \"This rise is approximately 1.5 times greater than the global average increase.\" I find 1.5 times greater to be somewhat ambiguous. I would propose to rephrase \"The rise is approximately 50% greater than the global average increase.\"\n* L59: \"focusing on extreme warming event\" -> plural events. At some places, the choice of plural versus singular forms should be checked before final publication, but that may be also something for the proofs.\n* L104: 'yet remain evident'\n* L111ff: \"the pattern of low-oxygen extreme events broadly resembles that of heatwave events, with a pattern correlation coefficient over 0.8 (p < 0.01).\" Maybe it would be good to be more specific here, writing something like 'the pattern of low-oxygen extreme event frequency and duration broadly resembles that of marine heatwave events, with pattern correlation coefficients over 0.8' (I presume it is two correlation coefficients)\n* L132: 'at each grid' do you mean something like 'in each grid cell'?\n* Fig S4: maybe you want to add to the caption of Fig S4 that panel b is the same as Fig 2a in the main text, may facilitate building bridges to the main text.\n* L165f: \"show a considerable rise in the annual CLH days across more than 70% of the world’s oceans, with an average increase exceeding 20 days\" do the 70% refer to the fraction of the ocean that experiences any increase in CHL events? proposal, sth like: \"show a considerable rise in annual CLH days, with increases in more than 70% of the world's ocean and a global average increase exceeding 20 days.\"\n* L171: something to double check: how can the global increase over 15 years be larger than 20 days (L165f), while the linear trend is about 1 per year globally, which would propose 15 days increase in 15 years.\n* L188 'the strengthening'\n* L190-192: not sure whether I understand 'a majority role in timeline of low-oxygen extreme events', maybe rephrase\n* L217-220: I am not sure I understand the sentence starting with 'This asymmetry..'. I would propose to rewrite the sentence. The finding of increased fish biomass loss where also CLH event days increase more is very interesting. It is clear that no causal relationship \"CLH event increase -> decline in fish biomass\" is necessarily implied, but it may open interesting scientific questions for future research.\n* L283f: to clarify \"However, the influence of NAO is relatively minor compared to other modes; hence it is not drawn in the figures.\" NAO is shown on the map in panel c, it's just not represented by a separate bar plot subpanel, right?\n* L284ff: I struggle to match this sentence with Fig. S8a, there does not seem to be an influence close to Greenland. Please double-check\n* L297: \"though we acknowledge that some signals (such as aerosols and volcanic activities) are still mixed in\" lacks some clarity, are you expressing they influence both anthropogenic and internal variability components? maybe rephrase\n* L334: Should be \"a marked increase in\", I think.\n* L433-434: \"Nevertheless, the simulated distribution is in good agreement with it in observations from WOA 2018\" WOA or WOD? in the beginning of the section you refer to WOD. To me it is not clear whether you compare the model against WOD or WOA (which is of course mainly building on WOD) in this section, or compare to both as indicated by the section title?"}, "author_response": "Response to Reviewer #1’s Comments\nI have re-read the manuscript and assessed the changes (as discussed in the response\nto reviewers). The authors have responded to my second review comments. I found\nthe manuscript in very good shape and I am satisfied with this revised version. The\ncaptions to each of the figures are clear. The distinction between trends and\nattribution to climate change is now accurate and clear. The manuscript reads very\nwell, no further changes from me.\nResponse: Thank you for your thorough re-evaluation of our manuscript and for your\npositive feedback on the revisions we have made. We are pleased to hear that you found\nthe manuscript in very good shape and overall readability have met your expectations.\nWe greatly appreciate your constructive comments throughout the review process,\nwhich have been instrumental in improving the quality and clarity of our work.\nResponse to Reviewer #2’s Comments\nThe authors have very carefully responded to all my comments. Thank you! I have\none final science related question, probably beyond the scope of the article: Why is\nthe rise in CHL events larger in high-biomass regions? It would be very interesting\nto provide some interpretation explaining this finding.\nResponse: Thank you for your positive feedback and for raising an interesting scientific\nquestion. We think that the larger rise in CLH events in these regions can be attributed\nto several factors. High-biomass areas are primarily located in near-coastal regions.\nThese areas typically have more abundant nutrient availability, which can support more\nfrequent and intense phytoplankton blooms, including those associated with CLH\nevents. Moreover, Guo et al. (2022) revealed that due to boundary currents, coastal\nprocesses, and ocean eddy fluxes, the heatwaves experienced in high-biomass regions\nmay be more severe than expected, which also provides a possible explanation for your\nquestions.\nWe agree that providing a detailed interpretation of this finding would be very\ninteresting and could be a valuable direction for future research. Although this is beyond\nthe scope of our current manuscript, we appreciate the suggestion and will consider it\nfor further investigation.\nThank you again for your insightful comments and contributions to our work.\nReference\nGuo, X., Gao, Y., Zhang, S., Wu, L., Chang, P., Cai, W., et al. (2022). Threat by marine\nheatwaves to adaptive large marine ecosystems in an eddy-resolving model. Nature\nClimate Change, 12(2), 179–186.\nI support publication of the manuscript as is. However, reading it once more, I have\na couple of small text-related comments:\nResponse: Thank you for your support for the publication of our manuscript. We\nappreciate your attention to detail and are grateful for your feedback. We have carefully\nrevised our manuscript and addressed the comments you have provided below to\nimprove the clarity and quality of the manuscript.\n* L209: \"This rise is approximately 1.5 times greater than the global average\nincrease.\" I find 1.5 times greater to be somewhat ambiguous. I would propose to\nrephrase \"The rise is approximately 50% greater than the global average increase.\"\nResponse: The sentence has been revised to: \"The rise is approximately 50% greater\nthan the global average increase.\"\n* L59: \"focusing on extreme warming event\" -> plural events. At some places, the\nchoice of plural versus singular forms should be checked before final publication,\nbut that may be also something for the proofs.\nResponse: Done.\n* L104: 'yet remain evident'\nResponse: Done.\n* L111ff: \"the pattern of low-oxygen extreme events broadly resembles that of\nheatwave events, with a pattern correlation coefficient over 0.8 (p < 0.01).\" Maybe it\nwould be good to be more specific here, writing something like 'the pattern of low-\noxygen extreme event frequency and duration broadly resembles that of marine\nheatwave events, with pattern correlation coefficients over 0.8' (I presume it is two\ncorrelation coefficients)\nResponse: We have revised this sentence following your suggestion.\n* L132: 'at each grid' do you mean something like 'in each grid cell'?\nResponse: Yes, it should be ‘in each grid cell’, we have corrected it in the revised\nmanuscript.\n* Fig S4: maybe you want to add to the caption of Fig S4 that panel b is the same as\nFig 2a in the main text, may facilitate building bridges to the main text.\nResponse: Thank you for your valuable suggestion. We have added this information to\nthe caption of Fig. S4.\n* L165f: \"show a considerable rise in the annual CLH days across more than 70%\nof the world’s oceans, with an average increase exceeding 20 days\" do the 70%\nrefer to the fraction of the ocean that experiences any increase in CLH events?\nproposal, sth like: \"show a considerable rise in annual CLH days, with increases in\nmore than 70% of the world's ocean and a global average increase exceeding 20 days.\"\nResponse: Yes, the 70% refer to the fraction of the ocean that experiences any increase\nin CLH events. Following your suggestion, this sentence has been revised to “show a\nconsiderable rise in annual CLH days, with increases in more than 70% of the world's\nocean and a global average increase exceeding 20 days.”\n* L171: something to double check: how can the global increase over 15 years be\nlarger than 20 days (L165f), while the linear trend is about 1 per year globally, which\nwould propose 15 days increase in 15 years.\nResponse: We have carefully checked the statement here. In line 165, the change is\ncalculated based on the difference between the average of two 15-yr period 1985-1999\nand 2000-2014 (that is, avg minus avg ), which is therefore exceeding\n_2000-2014 _1985-1999\n20 days. This does not contradict the result (linear trend of 0.96 day per year) described\nin L171.\n* L188 'the strengthening'\nResponse: Done.\n* L190-192: not sure whether I understand 'a majority role in timeline of low-oxygen\nextreme events', maybe rephrase\nResponse: We apologize for any confusion caused by the original phrasing. To improve\nclarity, we have rephrased the sentence as follows: “The joint occurrences (i.e., CLH\nevents) has become the dominant concern in the evolution of low-oxygen extreme\nevents”\n* L217-220: I am not sure I understand the sentence starting with 'This asymmetry..'.\nI would propose to rewrite the sentence. The finding of increased fish biomass loss\nwhere also CLH event days increase more is very interesting. It is clear that no causal\nrelationship \"CLH event increase -> decline in fish biomass\" is necessarily implied,\nbut it may open interesting scientific questions for future research.\nResponse: The term “asymmetry” here refers to the greater increase of extreme events\nin high-biomass areas compared to low-biomass areas, illustrating an asymmetrical\nincrease pattern between these regions. To enhance clarity, we have rewritten the\nsentence as follows: “This greater increase of CLH events in high-biomass regions\nimplies more severe damage to marine life and ecosystem”. Our discussions about why\nthe rise in CLH events is larger in high-biomass regions can be found in our response\nto Comment #1.\n* L283f: to clarify \"However, the influence of NAO is relatively minor compared to\nother modes; hence it is not drawn in the figures.\" NAO is shown on the map in panel\nc, it's just not represented by a separate bar plot subpanel, right?\nResponse: Yes, you are totally right.\n* L284ff: I struggle to match this sentence with Fig. S8a, there does not seem to be\nan influence close to Greenland. Please double-check\nResponse: We have carefully checked this sentence. The NAO indeed has some impact\non the oceans east of Greenland, as shown in Fig. R1 (despite the affected area is\nrelatively limited). Nevertheless, to avoid any confusion, we have removed this\nsentence.\nFigure R1. Percentage change in annual CLH event days in NAO positive phase. This\nfigure is same as Fig. S8a in the Supplementary material.\n* L297: \"though we acknowledge that some signals (such as aerosols and volcanic\nactivities) are still mixed in\" lacks some clarity, are you expressing they influence\nboth anthropogenic and internal variability components? maybe rephrase\nResponse: We apologize for the lack of clarity in this sentence. Here we intend to\nindicate that our method exclusively separates extreme events into anthropogenic and\ninternal variability components; it does not identify specific signals such as aerosols\nand volcanic activities. We have rephrased this sentence as follows: “our method does\nnot identify specific signals such as aerosols and volcanic activities”.\n* L334: Should be \"a marked increase in\", I think.\nResponse: Thank you for your comment. We have revised this sentence.\n* L433-434: \"Nevertheless, the simulated distribution is in good agreement with it in\nobservations from WOA 2018\" WOA or WOD? in the beginning of the section you\nrefer to WOD. To me it is not clear whether you compare the model against WOD or\nWOA (which is of course mainly building on WOD) in this section, or compare to\nboth as indicated by the section title?\nResponse: As described in Line 343 and 362, both WOD and WOA are used in model\nvalidation. We compare the model results against to both of them as indicated by the\nsection title."}]

Decision letter