Human-induced intensified seasonal cycle of sea surface temperature

Domain: earth_science
Published: 2024-05-10
Editorial outcome: major_revision

Abstract
Changes in the seasonal cycle of sea surface temperature (SST) have far-reaching ecological and societal implications. Previous studies have found an intensified SST seasonal cycle under global warming, but whether such changes have emerged in historical records remains largely unknown. Here, we reveal that the SST seasonal cycle globally has intensified by 3.9 ± 1.6% in recent four decades (1983–2022), with hotspot regions such as the northern subpolar gyres experiencing an intensification of up to 10%. Increased greenhouse gases are the primary driver of this intensification, and decreased anthropogenic aerosols also contribute. These changes in anthropogenic emissions lead to shallower mixed layer depths, reducing the thermal inertia of upper ocean and enhancing the seasonality of SST. In addition, the direct impacts of increased ocean heat uptake and suppressed seasonal amplitude of surface heat flux also contribute in the North Pacific and North Atlantic. The temperature seasonal cycle is intensified not only at the ocean surface, but throughout the mixed layer. The ramifications of this intensified SST seasonal cycle extend to the seasonal variation in upper-ocean oxygenation, a critical factor for most ocean ecosystems.

Introduction
Sea surface temperature (SST) plays a pivotal role in linking the atmosphere and ocean, attracting extensive attention to deciphering variations in annual and seasonal mean SST1–3. However, even in the absence of annual mean state changes, SST retains the potential to exert substantial climatic influence by modulating the amplitude of its seasonal cycle. These changes in the SST seasonal cycle can have profound climatic impacts, influencing phenomena such as marine heatwaves, Asian monsoons, precipitation, and the El Niño-Southern Oscillation4–7. Furthermore, changes in SST seasonality bear considerable ecological implications, particularly concerning oceanic oxygen content, a factor directly influencing ocean productivity and biogeochemical cycles8–12.

In a warmer future climate, climate models consistently project a global intensification of the SST seasonal cycle13–17. Under extreme warming scenarios, Coupled Model Intercomparison Project phase 5 (CMIP5) models indicate that the amplitude of the SST seasonal cycle could intensify by ~30% by the end of the 21st century compared to the pre-industrial levels14. The precise mechanisms propelling this intensification remain a subject of debate. Some studies have emphasized the thermodynamic effects of atmospheric circulation changes, which modify the seasonal cycle of surface heat fluxes and subsequently impact the seasonal cycle of SST18–20. Other studies have highlighted the crucial role of changes in the mixed layer depth (MLD) in driving the intensified SST amplitude13,14,17. As the Earth warms, increased ocean heat uptake leads to enhanced stratification in the upper ocean, resulting in a shallower MLD21–23. The reduced heat content within the shallower MLD diminishes thermal inertia, intensifying the SST’s seasonal cycle.

Despite the consensus on the intensified SST seasonal cycle in future warming scenarios, it remains challenging to determine whether such changes have already manifested themselves in the observations. Unlike future projections, the sparse observational record and the influence of internal variability make it much more difficult to detect these changes in the observational period. In addition, existing studies that have explored the detection of seasonal cycle changes in surface temperature have focused primarily on land24–28, leaving the emergence of surface temperature trends over the oceans largely elusive. Our study fills this gap by using a combination of observations and model simulations, unveiling a substantial enhancement of the SST seasonal cycle over most of the oceans during the past four decades.

Results
Intensified global SST seasonal cycle in observations
The amplitude of the SST seasonal cycle is assessed by calculating the difference between the maximum and minimum values of the annual cycle at each grid point (see “Amplitude of seasonal cycle” in Methods). For our investigation, we rely primarily on three widely-used observational SST datasets: the Extended Reconstructed SST dataset (ERSST) version 5, the Hadley Centre Sea Ice and SST dataset (HadISST) version 1.1, and the high-resolution Optimum Interpolation Sea Surface Temperature (OISST) version 2. These datasets cover the period from 1982 to 2022, enabling a detailed examination of the SST seasonal cycle over the past four decades (see “Observational and reanalysis datasets” in Methods). Pre-1982 observations are excluded from our analysis due to the sparse observations, particularly during the deep winter season in high latitudes29, which may not adequately capture the SST seasonal cycle. We denote the mean of the three datasets as OBS.

Climatologically, the Northern Hemisphere SST in OBS exhibits a maximum in August and a minimum in February (Fig. 1a), while the Southern Hemisphere shows the opposite evolution (Supplementary Fig. 1). Over the past four decades, the main changes in SST seasonality are an enhanced difference between the annual maximum and minimum temperature (Fig. 1a, bars), demonstrating a substantial increase in amplitude. The intensified SST seasonal cycle is further confirmed by the time series of the seasonal cycle amplitude in all three observations (Fig. 1c, black lines). Overall, the intensification in OBS corresponds to a statistically significant increase of 0.16±0.07°C (the confidence interval is 5–95%) from 1983 to 2022, representing a relative increase of 3.9±1.6% compared to the climatological mean of 4.2°C during 1983–1992 (Fig. 1d, gray bars). Note that the first year, 1982, is excluded from the trend analysis due to its extremely low value, which would skew the calculated trend.

To tease out the influence of internal variability and to derive an estimate of the forced response of the SST seasonal cycle to external forcing, we average the results of five simulation ensembles, comprising a total of 181 individual realizations (denoted as “EXT”; see “Climate model simulations” in Methods). The EXT ensemble successfully reproduces the observed seasonal cycle (Fig. 1b) and its amplification (Fig. 1c, red line), with an increase of ~3.4% over the period from 1983 to 2022. The intensification trend is significantly robust across all five simulation ensembles (Fig. 1d, red bars), ranging from ~1.7% in MPI-ESM1-LE to ~5.3% in CanESM5-LE. Moreover, this intensification trend appears in 35 out of 36 CMIP6 models, with intensifications ranging from 0.3% in GISS-E2-2-G to 7.0% in IPSL-CM6A-LR (Supplementary Fig. 2). The high consistency of global mean amplitude trends between observations and model simulations indicates that the intensification is primarily a robust climate response to external forcing, rather than being internally generated.

The linear trend of the observed amplitude between 1983 and 2022 displays considerable regional variations (Fig. 1e). The most substantial intensification occurs in the northern subpolar gyres, reaching ~11.9% between 45°N and 60°N. In the North Pacific, a distinctive horseshoe-shaped pattern of seasonal cycle enhancement emerges, while in the North Atlantic, the intensification primarily occurs in the western basin. In the Southern Hemisphere, prominent and extensive intensification (~5.0% between 35oS and 50oS) is observed north of the Atlantic Circumpolar Circulation (ACC). Intensifications also appear in the eastern tropical Pacific and northeast tropical Atlantic. However, specific regions exhibit a significant reduction in amplitude, including extensive areas of the western Pacific Ocean and the Indian Ocean, as well as the southern boundaries of the ACC. This leads to minor changes in the zonal mean amplitude near 30oN and 30oS. The pattern of amplitude change is consistent across the ERSST, HadISST, and OISST datasets (Supplementary Fig. 3).

By prescribing the external forcings, the EXT can well capture the major spatial features of observed changes in seasonal cycle amplitude, with a pattern correlation of 0.56. Positive trends are present in the North Pacific, North Atlantic, and the northern flank of the ACC, while muted and negative trends occur in the southern flank of the ACC, the western subtropical Pacific (comparing Fig. 1f, e). When considering only statistically significant trends in OBS (stippled regions in Fig. 1e), the pattern correlation increases to 0.85. This spatial resemblance between OBS and simulations remains consistent across all five simulation ensembles (Supplementary Fig. 4). Notably, discrepancies between observed and simulated amplitude changes are more pronounced in the Southern Hemisphere, particularly in the Pacific sector of the Southern Ocean. The pattern correlation is substantially lower in the Southern Hemisphere (0°−60°S) at 0.24 compared to 0.72 in the Northern Hemisphere (0°−60°N), suggesting the influence of large observational uncertainties in the data-sparse Southern Ocean. Furthermore, the low-frequency internal variability within the climate system30 may also contribute, as most of the discrepancy between EXT and OBS falls within the range of CMIP6 ensemble standard deviation (Supplementary Fig. 4b). Despite these discrepancies, the model simulations effectively replicate the temporal evolution and the magnitude of intensification trends of SST seasonal cycle in hotspot regions, thus highlighting the dominant role of external forcing.

Attribution of the observed trend in the SST seasonal cycle
To discern the effects of different external forcings, we employ model simulations from the Detection and Attribution Model Intercomparison Project (DAMIP) (see “Climate models simulations” in Methods), which enables separating the contributions of four major radiative forcings: greenhouse gases (GHGs), anthropogenic aerosols (AERs), stratospheric ozone (StratO3), and natural forcings (NATs) associated with volcanic and solar variations.

During 1982–2020, the continuous increase in GHGs concentration has amplified the SST seasonal cycle by approximately ~2.0% (Fig. 2a, red line), accounting for ~55.6% of the total historical enhancement in CMIP6 models (Fig. 2a, black line). The AERs also play a role in amplifying SST seasonality by ~0.8% (Fig. 2a, light blue line), as a result of subsequent reduction in aerosol emissions in Europe, North America, and China since the 1980s31. In contrast, contributions from StratO3 have barely any contribution to the seasonal amplitude changes (Fig. 2a, yellow line). The ~0.4% intensification caused by NATs largely results from two abrupt drops occurring in 1983 and 1992 and thus is not statistically significant (Fig. 2a, purple line). These drops follow the explosive volcanic eruptions of El Chichon in 1982 and Mount Pinatubo in 1991, respectively, underscoring the discrete impact of NATs on specific periods of the SST seasonal cycle.

Similar to its impact on temporal evolution, GHGs are the primary driver shaping the spatial pattern of amplitude changes in the historical simulations (comparing Fig. 2b to Fig. 1f), including the widespread intensification in the northern subpolar gyres and the dipole-like trends in the Southern Ocean. However, there are some differences between historical simulations and GHGs simulations. In the historical simulations, the intensification in the North Pacific is stronger and more concentrated within the subpolar gyre, whereas GHGs-induced intensification is weaker and spreads more into lower latitudes. The inclusion of AERs offers a plausible explanation for the differences between GHGs and historical simulations. For example, AERs generate an amplitude dipole in the North Pacific centered at ~35oN (Fig. 2c), where the amplitude increases on the poleward side and decreases on the equatorward side. This dipole structure contributes to the meridional contrast in the North Pacific in historical simulations, strongly enhancing the intensification in the subpolar regions. Overall, it is evident that the intensification of the SST seasonal cycle is predominantly driven by human-induced GHG emissions, and AERs also play a significant role in northern subpolar oceans.

Climate models robustly capture the observed intensification of the SST seasonal cycle, adding confidence to their future projections (Supplementary Fig. 5). In the high-emissions Shared Socioeconomic Pathway scenario 5–8.5 (SSP5-8.5), the SST seasonal cycle exhibits an almost linear upward trend, resulting in a remarkable ~10.6% enhancement by 2100 compared to the present-day level (mean of 2015–2024). Even under lower emission scenarios like SSP2-4.5 and SSP3-7.0, the intensification can reach ~4.8% and ~7.7%, respectively. It is important to note that certain ocean areas experience more considerable intensification due to regional variations. For instance, in the North Pacific and North Atlantic, the amplitude will enhance by ~10.3% between 45oN and 60oN, even in the modest SSP2-4.5 scenario (Supplementary Fig. 5b).

Mechanisms of the intensified SST seasonal cycle
To explore the mechanisms driving changes in the SST seasonal cycle over the past four decades, we utilize CMIP6 simulations and conduct a mixed layer budget analysis (see “Mixed Layer Heat Budget” in Methods). This approach enabled us to quantitatively assess the contributions of different factors (Figs. 3, 4), including the thermal forcing term Qm that is proportional to the ratio between the net surface heat flux (SHF) and MLD, as well as horizontal advection and other residual processes.

During 1983–2022, Qm exhibits positive trends across most of the global ocean (Fig. 3b), particularly noteworthy in the northern subpolar gyres, where induced seasonal intensification reaches 1.8°C north of 30°N, corresponding to a 58.1% increase relative to the 1983–1992 mean. The residual term emerges as the primary contributor to seasonal intensification in the equatorial oceans and also plays a role in the intensification on the northern flank of the ACC (Fig. 3d), and it strongly suppresses the seasonal amplitude in the Northern Hemisphere, indicating the importance of vertical entrainment or diffusive processes in these regions. In the global mean (Fig. 3a), the amplification of the SST seasonal cycle is predominantly attributed to positive contributions from Qm (~11.2%; red), and this effect is partially counteracted by the residual term (~−6.0%; orange). The contribution from the horizontal advection is deemed negligible (blue).

This investigation naturally prompts an inquiry into the primary factors shaping Qm. These factors include decreased MLD (Fig. 5a), reduced seasonal amplitude of MLD (Supplementary Fig. 6a), increased ocean heat uptake (Supplementary Fig. 6b), and enhanced amplitude of SHF across most of the global oceans (Supplementary Fig. 6c). Through detailed decomposition (Fig. 4; see “Mixed Layer Heat Budget” in Methods), our analysis highlights the predominant role of shallower MLD over the past four decades, largely attributed to enhanced upper ocean stratification32–35. The impact of decreased MLD on SST seasonality is most prominent in the northern subpolar oceans, where the induced seasonal intensification exhibits a remarkable 158.6% increase relative to the 1983–1992 mean. In the Southern Hemisphere, although the effect of MLD reduction is relatively smaller, it still contributes to an intensification around 35°S. Globally, the decreased MLD entirely drives the evolution in Qm, amplifying the SST seasonal cycle by 29.1% (Fig. 4e).

The intensification of the SST seasonal cycle due to decreased MLD is largely offset by the seasonal amplitude change in MLD (Fig. 4b), the latter is characterized by a broad reduction in mid-to high-latitude oceans (Supplementary Fig. 6a). The underlying mechanisms are depicted schematically in Supplementary Fig. 7: considering a scenario where SHF oscillates annually around zero and MLD also follows an annual-repeating evolution, albeit with reduced seasonal amplitude and no long-term trend. This results in a Qm characterized by decreased maxima in summer and increased minima in winter. However, the rate of decrease during summer surpasses the rate of increase in winter, owing to shallower MLD in warmer seasons. Consequently, the diminished seasonal cycle of MLD substantially contributes to a reduction in the SST seasonal cycle amplitude, effectively mitigating the contribution from decreasing MLD.

Changes in the annual mean and seasonal amplitude of SHF contribute to regional intensifications of SST seasonality (Fig. 4c, d): the anomalous heat uptake over the North Atlantic and around the Kuroshio (Supplementary Fig. 6b) intensify the SST seasonal cycle there, while increased amplitude of SHF over the North Pacific and the northern flank of ACC (Supplementary Fig. 6c) enhances the amplitude of the SST seasonal cycle. However, these contributions are considerably smaller compared to those resulting from decreased MLD (Fig. 4e).

Another line of evidence favoring the dominance of a shoaling mixed layer is the distribution of amplitude changes across the mixed layer (Fig. 5b). The intensification of the temperature seasonal cycle extends from the surface to the bottom of the mixed layer in all zonal bands. Particularly noteworthy are the changes in the northern subpolar gyres at around 55°N, where the regions correspond closely to areas of the strongest MLD shoaling (Fig. 5b, lines), further confirming the important role of a shallower mixed layer. Another intriguing feature is the amplitude reduction beneath the mixed layer, likely attributed to the overall shallowing of the mixed layer that impedes the propagation of the seasonal cycle signal into deeper ocean layers.

The changes in annual-mean MLD and zonal-mean temperature seasonality detected in climate models are also evident in the Institute of Atmospheric Physics (IAP) (Fig. 5c, d) and Ishii datasets (Supplementary Fig. 8a, b). However, climate models exhibit a deeper-reaching intensification in the northern high latitudes compared to the observations, which is consistent with the deeper climatological MLD in climate models. Notable differences exist in another observational dataset EN4, particularly in the Southern Hemisphere (Supplementary Fig. 8c, d). This discrepancy may be attributed to EN4’s use of temperature and salinity relaxation to climatology in areas where data is scarce36. Despite this, the strong consistency among results from IAP, Ishii, and model simulations provides further confidence in the robustness of our findings.

Implications of the intensified SST seasonal cycle
There is a close relationship between oxygen solubility and SST, wherein the oxygen solubility decreases with rising temperature. Hence, there has been a significant ocean deoxygenation trend since the mid-20th century along with global warming10,37,38. However, the impact of global warming on ocean oxygen solubility may extend beyond annual mean trends. With the intensified SST seasonal cycle, the contrast of surface dissolved oxygen between the seasonal maximum in winter and minimum in summer increases (Fig. 6a), demonstrating an intensification at a rate of ~3.7% during 1983–2022 (Fig. 6b). Moreover, models with a more intensified SST seasonal cycle would also see a more intensified seasonal cycle of dissolved oxygen, demonstrating a correlation of 0.84 (Fig. 6c). The spatial patterns of the intensification in these two seasonal cycles (Fig. 1f versus Fig. 6d) are also similar, with a pattern correlation of 0.60. This collectively underscores the impact of the intensification of the SST seasonal cycle on ocean oxygen solubility. Regional analysis reveals a substantial 15.7% enhancement in the seasonal contrast of dissolved oxygen in the subpolar North Atlantic (40–60oN, 60oW-0). Intensifications in the subpolar North Pacific (40–60oN, 140oE-120oW) and Southern Ocean (35–50oS) also reach 5.5% and 6.9%, respectively. These enhancements may foster the occurrence of hypoxic conditions characterized by exceedingly low oxygen concentrations39 by superimposing upon the long-term oxygen depletion trend.

Apart from dissolved oxygen, the air-sea CO2 flux (FCO2) is significantly influenced by SST variations. Our analysis reveals a notable intensification of FCO2 seasonality during 1983–2022 (Supplementary Fig. 9a). However, its spatial pattern differs from that of intensified SST seasonality. For instance, the intensification in the FCO2 seasonal cycle is notable in subtropical oceans and the southern flank of the ACC, where changes in the SST seasonality are relatively minor or even negative. Conversely, in regions such as the subpolar North Pacific north of 45°N, which experience a remarkable SST seasonal cycle enhancement, the FCO2 seasonal cycle is suppressed. Furthermore, the intermodel correlation between global mean changes in FCO2 and SST seasonal cycles is weak (r=0.14; Supplementary Fig. 9b), suggesting that mechanisms responsible for the intensified FCO2 seasonal cycle involve other factors, such as surface winds40.

Methods
Observational and reanalysis datasets
In this study, we focus on the amplitude changes of the SST seasonal cycle from 1982 to 2022. We utilize three observational datasets: (1) Extended Reconstructed SST (ERSST) version 5, covering the period from 1854 to the present with a 2°×2° grid resolution45, (2) Hadley Centre Sea Ice and SST dataset (HadISST) version 1.1, covering the period from 1870 to the present with a 1°×1° grid resolution46, and (3) Optimum Interpolation Sea Surface Temperature (OISST) version 2, covering the period from 1982 to the present with a 0.25°×0.25° grid resolution47.

To examine changes in the seasonal cycle of mixed layer temperature and the annual mean MLD from 1982 to 2022, we use three monthly objective analysis ocean dataset: (1) the Institute of Atmospheric Physics (IAP) ocean temperature and salinity analysis, covering the period from 1940 to the present with a 1°×1° grid resolution and 41 vertical levels for the upper 2,000 m48, (2) the Ishii v7.3, covering the period from 1955 to 2019 with a 1°×1° grid resolution and 28 vertical levels for the upper 3000 m49, and (3) the EN4 data, covering the period from 1900 to present with a 1°×1° grid resolution and 42 vertical levels extending to beyond 5,000 m36. It is important to note that the EN4 dataset applies temperature and salinity relaxation to climatology in regions with sparse data coverage, and thus conducting long-term trend analysis with EN4 could potentially yield problematic results36. To ensure consistent analysis, all the above datasets were interpolated onto a 2°×2° grid using bilinear interpolation.

Climate model simulations
We use the CMIP6 model outputs of SST, ocean temperature, zonal velocity, meridional velocity, vertical velocity, surface heat flux, MLD, surface dissolved oxygen, surface CO2 flux from the historical and SSP scenarios simulations based on data availability (Supplementary Table 1)50. Since CMIP6 historical simulations ended in 2014, we extend them to 2022 using SSP5-8.5 scenarios51. Only the first member of each model is used to ensure equal weight in the multi-model ensemble mean analysis. In addition, we use historical and SSP585 SST data from four large ensemble simulations (LENS): (1) 40 ensemble members of the ACCESS-ESM1-552, (2) 25 ensemble members of the CanESM5-LE53, (3) 30 ensemble members of the MPI-ESM1-2-LR54, and (4) 50 ensemble members of the MIROC655. For the analysis of MLD and surface dissolved oxygen, only the first three LENS are utilized due to data availability.

To assess the contributions of different anthropogenic and natural forcings to the SST seasonal cycle changes, we use SST output from single forcing simulations from the Detection and Attribution Model Intercomparison Project (DAMIP)56. These simulations include scenarios with GHGs only, AERs only, StratO3 only, and NATs only simulations. We select seven models that have more than three members for single forcings simulations, and the simulations cover the period of 1982–2020 (Supplementary Table 1). All available members from each model are included in the analysis. All the outputs are interpolated onto a 2°×2° grid using a bilinear interpolation.

Amplitude of seasonal cycle
The amplitude of the SST seasonal cycle is determined at each grid point by calculating the difference between the monthly maximum and minimum SST values. Although the timings of maximum and minimum SST can change under climate change, the magnitude of this change is relatively small, typically only a few days. Therefore, we obtain the months of maximum and minimum SST based on the climatology spanning 1970–2000 from the ERSST. These identified months serve as the basis for calculating the amplitude of the seasonal cycle for each year in both observations and simulations. A similar method is applied to surface heat flux, ocean temperature, surface dissolved oxygen, and air-sea CO2 flux data to determine their respective seasonal cycle amplitudes. For ocean temperature and MLD, the climatology is based on the period 1970–2000 from the IAP dataset. For surface heat flux, surface dissolved oxygen and air-sea CO2 flux, the climatology is based on the same period from the CMIP6 multi-model mean.

Mixed layer depth
CMIP6 models provide the variable mlotst, which is the MLD calculated instantaneously on the model timestep. It is defined as the depth at which the potential density exceeds the sea surface density at a criterion of 0.125 kg m−3. In comparison, our calculations for the IAP dataset utilize monthly potential density with a 0.01 kg m−3 criterion57. The potential density is calculated from temperature and salinity fields using the TEOS-10 Gibbs Seawater (GSW) toolbox. While a commonly used criterion of 0.03 kg m−3 would yield similar results, the 0.01 kg m−3 criterion was chosen due to its perfect alignment with the depth at which the intensification in the temperature seasonal cycle can be observed (Fig. 5b, d). Our analysis of MLD based on the IAP, Ishii, and EXT simulations shows similar decreasing trend from 1983 to 2022. However, there is a notable discrepancy observed with the EN4 dataset, possibly due to the fact that EN4 incorporates temperature and salinity relaxation to the climatology in data-sparse regions36, while IAP and Ishii do not employ such a relaxation method.

Mixed layer heat budget
We investigated the heat budget of the ocean mixed layer as the following:

∂Tm/∂t = Qm − um · ∇Tm + r, (1)

Here, Tm represents the mixed layer temperature. The thermal forcing term Qm is determined by the surface heat flux Qnet and the mixed layer depth h. The constants Cp and ρ are the density and heat capacity of seawater, respectively. um · ∇Tm represents the horizontal advection term, where um is the horizontal velocity averaged over the mixed layer. The vertical entrainment and other unresolved processes are calculated as a residual term r. We then evaluated the contributions of different factors to changes in Tm by integrating each tendency term on the right-hand side of the equation over time from January 1982 to December 2022. To obtain the amplitude of the time-integrated tendency terms, we use values corresponding to the months of maximum and minimum Tm.

In addition, we decomposed the thermal forcing term into contributions from four factors14: annual mean changes in Qnet, seasonal cycle amplitude changes in Qnet, annual mean changes in MLD, and seasonal cycle amplitude changes in MLD. To assess the impact of each factor, we calculated the time-integrated thermal forcing term with one time-evolving factor at a time, while keeping the other three factors constant at their climatological means during the integration. Due to data availability (Supplementary Table 1), the heat budget analysis involves eight models: CAMS-CSM1-0, CAS-ESM2-0, CanESM5, CESM2-WACCM, FGOALS-f3-L, IPSL-CM6A-LR, NESM3.

Statistical significance test. The statistical significance of the correlation coefficients is determined by using an “effective sample size” N*58:

N* = N (1 − r1 r2) / (1 + r1 r2), (2)

Where N is the number of available time steps and r1 and r2 are lag-one autocorrelation coefficients of each variable.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "General Comments:\nThis study examined the intensificafion of SST seasonal cycle using observafional data and CMIP6/DAMIP\nmodels output. The authors suggested 4.6 (±1.5)% intensificafion of SST seasonal cycle’s amplitude\nduring 1960-2022, especially, northern subpolar gyre and northern flank of the ACC, where the\namplificafion is up to 10%. CMIP models show similar paftern of SC amplificafion and employing DAMIP\noutputs, anthropogenic GHG is more responsible for this reinforcement than aerosol and natural climate\nvariability. Addifionally, they show indicafion of its impacts on ocean biogeochemical properfies like\ndissolved oxygen and ecosystem.\nIn general, the manuscript is well wriften and reads quite well. I found this study quite interesfing and\npotenfial to share more widely with interdisciplinary fields. On the other hand, I have some quesfions\nand suggesfions for adding more discussions. After the authors consider my points, I would recommend\nthis work published in Nature Communicafions.\nSpecific Comments.\n1. Figure 3d and f, zonal-mean of MLD from observafion and models. Throughout the manuscript, the\nauthors emphasize the MLD shoaling in the higher lafitude in NH and SH (subpolar gyre and ACC). I could\nsee some remarkable distribufion of shoaling from 2D plots. However, the lines in Figs.3d and f look\nquite different between observafion and models. In the NH (let say, 50-60N), observafion (IAP) shows\nquite a small difference between two periods while model (EXT) shows a clear difference. In the SH, the\ndifference between two periods does not look similar between IAP and EXT. Also, climatological MLD of\nIAP seems too shallow or very different from EXT (40m and 100m in NH and 60m and 100m). I am\nwondering how reliable IAP is as a benchmark. And given that IAP is a qualified reference for models, I\nam curious if the climatological bias of EXT’s MLD influences interpretafion of results. I would\nrecommend that some argument of this point can be added.\n2. While I see a good agreement between SC intensificafion and MLD shoaling especially, subpolar\nregion. However, I am curious about the mechanism of MLD shoaling and its contribufion to SC\nintensificafion in those regions. The authors aftribute that to the GHG and corresponding ocean heat\nuptake using DAMIP and referring literatures, but global warming influence dynamics of atmosphere and\nocean as well. To me, the authors’ argument sounds like one dimensional mechanism. However, in\nparficular, subpolar gyre and ACC are more complicated because of ocean current and sea ice formafion\n(also maybe atmospheric circulafion). I am wondering if the changes in those components can also\nchange MLD and seasonal cycle of SST. Actually, the authors menfion ACC change due to atmospheric\ncirculafion, but how about in the NH, Pacific and Atlanfic (for example, Gulf Stream and\nKuroshio/Oyashio changes and influencing SST SC and MLD)?\n3. Biogeochemical aspect. While the authors show implicafion for dissolved oxygen and marine\necosystem, I am curious how air-sea carbon dioxide flux would respond to SC intensificafion. Because\nsubpolar region of NH and SH is a sink of air-sea CO2 flux, it is important in global carbon cycle. I totally\nunderstand there is limitafion of space, but really interesfing to see any discussion or implicafion about\nimpacts on air-sea CO2 flux. The authors would NOT need to add any plots for this point. I would like to\nsee how the authors think about this.\n4. Personally, the suppression of SST SC in 1983 and 1992 due to volcanic erupfion is quite interesfing.\nI would suppose that this suppression can be aftributed to less downwelling shortwave radiafion because\nof abundant aerosol emission. Would the authors see any agreement of MLD deepening and shortwave\nradiafion change in 2D distribufion? On the other hand, GHG contribufion to intensificafion of SST SC is\ndue to more re-emission of longwave radiafion down to the surface. So, I am wondering if the\ndownwelling longwave radiafion at the surface and MLD shoaling match well. Perhaps, this might\nstrengthen the authors’ discussion on causality of SST SC intensificafion due to MLD shoaling.\nMinor Comments.\n1. Fig1a and b. The capfion is not clear.\n2. Line 133. Blue => light blue?\n3. Fig.2a. The thickness of line in the legend can be widen a bit more?\n4. Line 136. Europe, North America, and China?\n5. Lines 149-150. But, not significant.", "reviewer_2": "Review of Human-induced intensified seasonal cycle of sea surface temperature by Liu et al.\nIn this study, the authors use several observafional data products and CMIP6 model results to invesfigate\nthe intensificafion of the SST seasonal cycle from 1960 to the present. The results show that there has\nbeen a 4.6% increase in the SST seasonal cycle in the global ocean, with the magnitude of the increase\nvarying widely among ocean regions. Using climate model output from DAMIP, they conclude that most\nof the increase in SST seasonal cycle amplitude has been driven by anthropogenic greenhouse gas\nemissions. The authors argue that the increase in SST seasonal cycle amplitude is mainly due to a\ndecrease in the annual mean MLD. The authors also suggest that the increase in the SST seasonal cycle\nmay have affected the seasonal cycle of sea surface dissolved oxygen through changes in its surface\nsolubility due to the amplificafion of the SST seasonal cycle.\nAs the authors menfioned in the introducfion, there have already been many studies that have proposed\nthe concept of long-term changes in the SST seasonal cycle itself using observafional data, as well as\nstudies that have used climate models to predict future changes in the amplitude of the SST seasonal\ncycle. Therefore, the novelty of this study lies only in the spafial coverage of the assessment of historical\nlong-term trends. The method used for aftribufing changes in the SST seasonal cycle is an established\none (using exisfing outputs from DAMIP). The authors' conclusion on the mechanism is also not a new\nhypothesis, as it has been proposed in previous studies on the predicfion of future changes in SST\nseasonal cycle amplitudes using climate models. The implicafions for changes in the seasonal amplitude\nof sea surface dissolved oxygen may be a new aspect of this study, but the validity of the model\ndissolved oxygen data used has not been fully invesfigated.\nThroughout the manuscript, I have three main concerns and several more to follow.\nMajor concerns:\n(1) Reliability of SST data before 1982\nAs the authors know, the number of SST observafions has changed dramafically since the advent of the\nsatellite observafions in 1982, i.e., the number of pre-1982 observafions is remarkably small. It is a well-\nknown fact that the number of observafions is fatally low, especially in the Southern Ocean and in the\ndeep winter season in high lafitudes, and that there are not enough observafions to fully resolve the SST\nseasonal cycle (see, for example, the spafial coverage and number of observafions archived in the World\nOcean Database). I have doubts about the inclusion of the pre-1982 period in the analysis of this study\nbased on some results read from the figures in the manuscript, such as\n-Figure S1c-f: The winter trends in each hemisphere, which are important for defining the amplitudes of\nthe seasonal cycle, differ between products.\n-Figure S3: The magnitude and sign of the trend in the Southern Ocean varies greatly when the fime\nperiod is changed (1960- or 1982-) for both ERSST and HadISST (as also the authors menfion in L103 but\nassume as minor). Note that the consistency between the two products does not resolve this concern, as\nthe two products use similar data sources prior to 1982, just with different mapping methods.\n-Figure S1 a and b: There appears to be a large discrepancy between the model and observafions in\nFigure S1 a and b (for a fair comparison, please use the same y-axis range for these panels). This seems\nto be due to the strong quesfionable trend in the Indian Ocean sector of the Southern Ocean that\nappears in the calculafion since 1960. It is true that there could be a bias in the model, but if the\nobservafions contain errors, wouldn't it be befter to start from 1982 so that the model and observafions\nare more consistent?\n-Figure 1c: There seems to be no significant trend before 1980.\n(2) Conclusion about the mechanism\nThe authors concluded that the decrease in the annual mean MLD is the main \"mechanism\" for the\nincrease in the amplitude of the SST seasonal cycle by showing only an intermodel \"correlafion\"\nbetween changes in the annual mean MLD and SST amplificafion. However, this is only for the global\nmean trend, and the spafial correspondence of the signals is not very good, although the authors say it is\ngood (comparing Fig. 3c to 1e). Even in the zonal mean (Fig. 3d), there is no dominant mixed layer\nthinning north of 50N, where the seasonal amplificafion is greatest, and the seasonal SST amplitude is\nnot enhanced despite a strong mixed layer shallowing trend at 30S. Again, this may be due to errors in\nthe data product (in both IAP and Ishii) due to the lack of verfical profiles prior to the Argo era, as the\nmodel appears to be befter (Fig. 3f).\nSince there may be several other mechanisms for the seasonal amplificafion of the SST besides the one\nrejected by the authors (changes in the seasonal amplitude of the sea surface flux) (e.g. changes in the\nannual mean sea surface heat flux and changes in the seasonal mean MLD), I believe that the authors\ncan only point out the relafionship between the annual mean MLD and the SST amplificafion, and not\ndetermine the mechanism (as the authors also menfion in L203).\n3) Possibilifies for modeled dissolved oxygen\nIn general, even the Earth system models (ESMs) in the CMIP6 generafion suffer from serious dissolved\noxygen biases mainly due to incomplete parameterizafions of air-sea exchange process and biological\nprocesses. Although it is difficult to obtain the observafional trend of monthly sea surface dissolved\noxygen due to lack of historical observafional data, the model climatological seasonal cycle should be\nvalidated using recent climatological observafional data products, for example, World Ocean Atlas\n(hftps://www.ncei.noaa.gov/products/world-ocean-atlas) and Gridded Ocean Biogeochemistry from\nArfificial Intelligence – Oxygen (hftps://www.ncei.noaa.gov/access/metadata/landing-\npage/bin/iso?id=gov.noaa.nodc:0259304).\nOther important concerns:\nL34: In references 8-11, I could not find any reference to “the effect of changes in the SST seasonal cycle\non marine ecosystems\". Please elaborate this sentence.\nL82：What does the plus/minus range mean? Please clarify the stafisfical methodology in the text.\nL93: This is true only for the case of the Northern hemisphere (Fig. 1a and b). It seems that there are\nsome differences in the case of the Southern hemisphere (Fig. S3a and b)\nL113: How did the authors esfimate the degree of freedom for the stafisfical test? Using the number of\ndata gird point as the degree of freedom is overesfimate because adjacent grids are not necessarily\nindependent. Please clarify the stafisfical method.\nL116 “Some discrepancies between observafions and models may arise either from the influence of low-\nfrequency internal variability”: Could you esfimate the range of internal variability (ensemble standard\ndeviafion) using the ensemble members and compare the range with the discrepancies between\nobservafions and models?\nL215 “Another intriguing feature is the widespread amplitude reducfion beneath the mixed layer,\naftributed to the overall shallowing of the mixed layer that impedes the propagafion of the seasonal\ncycle signal into deeper ocean layers.”: Interesfing.\nL242: “deoxygenafion” rather than “oxygenafion”?\nL243 “The enhanced seasonal contrast of upper ocean oxygen level may foster the occurrence of hypoxic\ncondifions characterized by exceedingly low oxygen concentrafions by superimposing upon the long-\nterm oxygen deplefion trend over the past 60 years”: To quanfify this impact, please show the spafial\nmap of the relafive changes in the sea surface oxygen amplitude (i.e., oxygen amplitude changes over\nthe past 60 years divided by climatological seasonal amplitude).\nL325：”This calculafion is based on the climatology of the period 1970-2000 from the ERSST”: I could\nnot get this. Please elaborate this.\nL338: “The MLD derived with the 0.01 kg m-3 threshold aligns perfectly with regions that exhibit an\nintensified seasonal cycle of upper ocean temperature at all lafitudes”: I think this is result and thus this\nis not placed in the “Method”.\nFigure 1a and b: The left y-axis should be “Monthly trend” rather than “Amplitude trend”?\nFigure S6 and Method: Why did not the authors use observafional heat flux product (OAFlux\nhftps://climatedataguide.ucar.edu/climate-data/oaflux-objecfively-analyzed-air-sea-fluxes-global-oceans,\nand JOFURO3 hftps://www.j-ofuro.com/en/)? Observafional data is of limited durafion, but should be\nused for validafion of atmospheric reanalysis.\nFigure S6b: This analysis also should be done for each hemisphere like Fig. S7.\nL176: What about the spafial correspondence between the heat flux changes and SST amplitude\nchange?\nData: Is there any reason why the authors did not use the CESM2 large ensemble\n(hftps://esd.copernicus.org/arficles/12/1393/2021/)? The CESM2 is one of CMIP6 generafion models\nand it has outputs from single forcing experiment\n(hftps://journals.ametsoc.org/view/journals/clim/aop/JCLI-D-22-0666.1/JCLI-D-22-0666.1.xml).\nMethod: 0.01 kg m-3 is too small for the threshold when determining the MLD from the Grid data\nproduct. As noted in the reference, 0.01 kg m-3 is a value to use for raw CTD profiles, rather 0.1 kg m-3\nor 0.125 kg m-3 is more common because verfical profiles in gridded product are verfically smoothed.\nFigure S4a: It would be befter to show the spafial map of the CMIP6 model spread (standard\ndeviafion) next to Fig. S4a."}, "author_response": "We thank the reviewers for their insightful comments and detailed suggestions, which help\nimprove the manuscript. We have focused on satellite-era data (after 1982) to ensure the robustness\nof our results. In addition, we have conducted a comprehensive heat budget analysis to demonstrate\nthe primary role of decreased MLD in intensifying the SST seasonal cycle. We appreciate the\nvaluable input from both reviewers, which has motivated us to improve the quality of our work.\nIn the following, we provide a point-by-point response to both reviewers.\nNote that major changes are highlighted in red in the revised text.\nReviewer #1 (Remarks to the Author):\nGeneral Comments:\nThis study examined the intensification of SST seasonal cycle using observational data and\nCMIP6/DAMIP models output. The authors suggested 4.6 (±1.5)% intensification of SST seasonal\ncycle’s amplitude during 1960-2022, especially, northern subpolar gyre and northern flank of the\nACC, where the amplification is up to 10%. CMIP models show similar pattern of SC\namplification and employing DAMIP outputs, anthropogenic GHG is more responsible for this\nreinforcement than aerosol and natural climate variability. Additionally, they show indication of\nits impacts on ocean biogeochemical properties like dissolved oxygen and ecosystem.\nIn general, the manuscript is well written and reads quite well. I found this study quite interesting\nand potential to share more widely with interdisciplinary fields. On the other hand, I have some\nquestions and suggestions for adding more discussions. After the authors consider my points, I\nwould recommend this work published in Nature Communications.\nResponse: Thank you for your positive feedback and valuable suggestions. We have made\nsubstantial revisions to improve the manuscript. Specifically, we have performed a quantification\nanalysis using a mixed-layer heat budget to assess the relative contributions of various factors,\nwhich provides compelling evidence that decreased MLD is the primary driver of SST seasonality\nintensification. In addition, we have expanded our discussion to include an examination of the\ndiscrepancies between modeled and simulated MLD, and the implication of the intensified\nseasonal amplitude for surface CO2 flux.\nBelow are our point-by-point responses to your comments.\nSpecific Comments.\n1. Figure 3d and f, zonal-mean of MLD from observation and models. Throughout the manuscript,\nthe authors emphasize the MLD shoaling in the higher latitude in NH and SH (subpolar gyre and\nACC). I could see some remarkable distribution of shoaling from 2D plots. However, the lines in\nFigs.3d and f look quite different between observation and models. In the NH (let say, 50-60N),\nobservation (IAP) shows quite a small difference between two periods while model (EXT) shows\na clear difference. In the SH, the difference between two periods does not look similar between\nIAP and EXT. Also, climatological MLD of IAP seems too shallow or very different from EXT\n(40m and 100m in NH and 60m and 100m). I am wondering how reliable IAP is as a benchmark.\nAnd given that IAP is a qualified reference for models, I am curious if the climatological bias of\nEXT’s MLD influences interpretation of results. ·\nResponse: Thank you for bringing this to our attention. We have revised the analysis to present\ntrends in zonal-mean MLD over the period 1983-2022 (purple line in Fig. R1a), rather than the\nprevious average over 2013-2022. This adjustment helps to mitigate the influence of internal\nvariability on the results and provides a more accurate representation of long-term trends. As a\nresult, the discrepancy in MLD between the IAP observations and CMIP6 simulations in the 50-\n60N region has been significantly reduced, with both showing a notable decrease in MLD.\nHowever, the difference between observations and simulations in the SH is still large. One possible\nreason is the less reliable observational datasets in the SH. In addition, differences in the criteria\nused to calculate MLD between the IAP dataset and CMIP6 simulations may contribute to\ndiscrepancies. While CMIP6 calculates MLD instantaneously using a criterion of 0.125 kg/m3,\nour calculations for the IAP dataset utilize monthly temperature and salinity fields with a criterion\nof 0.01 kg/m3. This 0.01 kg/m3 criterion was chosen because the monthly fields are smoother\ncompared to the CMIP6 instantaneous fields, and the MLD defined in this way aligns well with\nthe intensified temperature seasonal cycle (blue line in Fig. R1b). If we were to use the 0.125\nkg/m3 criterion to define MLD in the IAP dataset, the results would indicate excessive depth,\nespecially at high latitudes (purple line in Fig. R1b), and would deviate significantly from the\nintensification trend.\nWe have incorporated this discussion on the potential causes for the discrepancy between observed\nand simulated MLD into the Method Section for clarity (Lines 354-361).\n2. While I see a good agreement between SC intensification and MLD shoaling especially,\nsubpolar region. However, I am curious about the mechanism of MLD shoaling and its\ncontribution to SC intensification in those regions. The authors attribute that to the GHG and\ncorresponding ocean heat uptake using DAMIP and referring literatures, but global warming\ninfluence dynamics of atmosphere and ocean as well. To me, the authors’ argument sounds like\none dimensional mechanism. However, in particular, subpolar gyre and ACC are more\ncomplicated because of ocean current and sea ice formation (also maybe atmospheric circulation).\nI am wondering if the changes in those components can also change MLD and seasonal cycle of\nSST. Actually, the authors mention ACC change due to atmospheric circulation, but how about in\nthe NH, Pacific and Atlantic (for example, Gulf Stream and Kuroshio/Oyashio changes and\ninfluencing SST SC and MLD)?\nResponse: We appreciate your perspective comment. To provide a comprehensive analysis of the\nunderlying mechanisms of the intensified SST seasonal cycle, we have conducted a mixed layer\nheat budget analysis. This approach allows us to quantitatively assess the contributions of various\nfactors, including the horizontal advection, the thermal forcing term, (influenced by net surface\nheat flux and MLD), and other residual processes (see Methods for more details).\nOur analysis reveals a limited role of ocean circulation changes in shaping the SST seasonality\nchanges, both in its temporal evolution (Fig. R2a, blue line) and trend pattern (Fig. R2c). In\naddition, the thermal forcing term is found to be dominant in driving the intensified SST seasonal\ncycle (Fig. R2a,b).\nFurthermore, we have decomposed the thermal forcing term into contributions from different\nfactors, including annual mean changes in SHF, seasonal cycle amplitude changes in SHF, annual\nmean changes in MLD, and seasonal cycle amplitude changes in MLD. This analysis reveals that\nthe decrease in MLD is the primary contributor (Fig. 4e), driving the overall intensification of the\nseasonal cycle.\nIn the revised manuscript, the section “Mechanisms of the intensified SST seasonal cycle” has\nbeen rewritten to include these findings (Lines 175-222), providing a compelling line of evidence\nfor the dominant role of decreased MLD.\n3. Biogeochemical aspect. While the authors show implication for dissolved oxygen and marine\necosystem, I am curious how air-sea carbon dioxide flux would respond to SC intensification.\nBecause subpolar region of NH and SH is a sink of air-sea CO2 flux, it is important in global\ncarbon cycle. I totally understand there is limitation of space, but really interesting to see any\ndiscussion or implication about impacts on air-sea CO2 flux. The authors would NOT need to add\nany plots for this point. I would like to see how the authors think about this.\nResponse: Thank you for your insightful comment. We have explored the response of air-sea CO2\nflux (FCO2) seasonal amplitude using the CMIP6 field fgco2 during the period 1982-2022.\nOur analysis reveals a significant intensification of FCO2 seasonality (Fig. R3a). However, its\nspatial pattern differs from that of the intensified SST seasonality. For instance, the intensification\nof the FCO2 seasonal cycle is notable in the subtropical oceans and the southern flank of the ACC,\nwhere changes in the SST seasonality are relatively small or even negative. Conversely, the FCO2\nseasonal cycle is suppressed in regions such as the subpolar North Pacific north of 45°N, which\nexperience a remarkable enhancement of the SST seasonal cycle. Furthermore, the intermodel\ncorrelation between global mean changes in FCO2 and SST seasonal cycles is low (r=0.1; Fig.\nR3b), indicating that the mechanisms responsible for the intensified FCO2 seasonal cycle involve\nother factors, such as surface winds and biological processes (Lerner et al. 2021).\nThe above discussion has been incorporated into the revised manuscript (Line 266-276&\nSupplementary Fig. 9) to provide a more comprehensive analysis of the biogeochemical\nimplications of the intensified SST seasonal cycle.\n4. Personally, the suppression of SST SC in 1983 and 1992 due to volcanic eruption is quite\ninteresting. I would suppose that this suppression can be attributed to less downwelling shortwave\nradiation because of abundant aerosol emission. Would the authors see any agreement of MLD\ndeepening and shortwave radiation change in 2D distribution? On the other hand, GHG\ncontribution to intensification of SST SC is due to more re-emission of longwave radiation down\nto the surface. So, I am wondering if the downwelling longwave radiation at the surface and MLD\nshoaling match well. Perhaps, this might strengthen the authors’ discussion on causality of SST\nSC intensification due to MLD shoaling.\nResponse: Thank you for raising this point. As you pointed out, changes in surface heat flux (SHF)\nand MLD, as well as changes in other ocean dynamical processes, may all contribute to changes\nin the SST seasonal cycle. To address this, we have conducted a comprehensive heat budget\nanalysis to quantitatively assess the contributions of different factors, as detailed in our response\nto your major comment #2. Our analysis reveals distinct mechanisms underlying the suppression\nof the SST seasonal cycle during the 1983 and 1992 volcanic events.\nSpecifically, we found that the suppression in 1992 is primarily induced by the thermal forcing\neffect alone, whereas the suppression in 1983 is influenced by both thermal forcing and ocean\ncirculation changes (Fig. R4a). Moreover, the controlling factors for the thermal forcing term in\nthe two events are also different: the 1992 event is mainly associated with changes in the MLD\n(Fig. R4c), while the 1983 event is largely contributed by amplitude changes in the SHF (Fig. R4b).\nGiven our focus on the long-term trend in SST seasonal cycle changes, we have deferred detailed\ninvestigations of these events to future analyses.\nRegarding the relationship between GHG-induced heat uptake and SST seasonal cycle\nintensification, our analysis indicates that although the ocean has absorbed more heat over the past\nfour decades, its direct impact on SST seasonal cycle intensification is limited (Fig. 4c, e). Instead,\nwe found that the decreased MLD plays a significant role in driving the intensification (Fig. 4a).\nIn the revised manuscript, we have thoroughly discussed the contribution of different factors to\nthe long-term intensification trend based on the mixed layer budget analysis, as outlined in Lines\n182-222.\nMinor Comments.\n1. Fig1a and b. The caption is not clear.\nResponse: Thank you, and we have updated the caption to clarify that the dotted lines represent\nthe SST seasonal cycle climatology, while the bars represent the trends.\n2. Line 133. Blue => light blue?\nResponse: Thank you, and it has been corrected.\n3. Fig.2a. The thickness of line in the legend can be widen a bit more?\nResponse: Done.\n4. Line 136. Europe, North America, and China?\nResponse: Revised as suggested.\n5. Lines 149-150. But, not significant.\nResponse: Thank you for pointing this out, and we have removed the description in the revision.\nReviewer #2 (Remarks to the Author):\nReview of Human-induced intensified seasonal cycle of sea surface temperature by Liu et al.\nIn this study, the authors use several observational data products and CMIP6 model results to\ninvestigate the intensification of the SST seasonal cycle from 1960 to the present. The results show\nthat there has been a 4.6% increase in the SST seasonal cycle in the global ocean, with the\nmagnitude of the increase varying widely among ocean regions. Using climate model output from\nDAMIP, they conclude that most of the increase in SST seasonal cycle amplitude has been driven\nby anthropogenic greenhouse gas emissions. The authors argue that the increase in SST seasonal\ncycle amplitude is mainly due to a decrease in the annual mean MLD. The authors also suggest\nthat the increase in the SST seasonal cycle may have affected the seasonal cycle of sea surface\ndissolved oxygen through changes in its surface solubility due to the amplification of the SST\nseasonal cycle.\nAs the authors mentioned in the introduction, there have already been many studies that have\nproposed the concept of long-term changes in the SST seasonal cycle itself using observational\ndata, as well as studies that have used climate models to predict future changes in the amplitude\nof the SST seasonal cycle. Therefore, the novelty of this study lies only in the spatial coverage of\nthe assessment of historical long-term trends. The method used for attributing changes in the SST\nseasonal cycle is an established one (using existing outputs from DAMIP). The authors' conclusion\non the mechanism is also not a new hypothesis, as it has been proposed in previous studies on the\nprediction of future changes in SST seasonal cycle amplitudes using climate models. The\nimplications for changes in the seasonal amplitude of sea surface dissolved oxygen may be a new\naspect of this study, but the validity of the model dissolved oxygen data used has not been fully\ninvestigated.\nThroughout the manuscript, I have three main concerns and several more to follow.\nResponse: We appreciate your perceptive comments, which have helped us improve the\nmanuscript substantially. As you can see from the revised version of our manuscript, considerable\nefforts have been made to address your comments, particularly in the following areas:\n(1) Use of data after 1982: To ensure the robustness of our results, we have focused on satellite-\nera data.\n(2) Quantification analysis based on a heat budget: We have conducted a mixed-layer heat\nbudget analysis to quantify the contribution of different factors to the intensification of the\nSST seasonal cycle. Through this analysis, we have identified decreased MLD as the primary\ndriver of SST seasonality intensification.\n(3) Validation of dissolved oxygen dataset: We have validated model simulations against the\nWOA18 to ensure reliability in simulating the climatological seasonal cycle of surface\ndissolved oxygen.\n(4) Internal variability: To assess the influence of internal variability on our results, we have\nestimated the model spread and found that most of the discrepancy falls within the range of\ninternal variability.\nBelow are our point-by-point responses to your comments.\nMajor concerns:\n(1) Reliability of SST data before 1982\nAs the authors know, the number of SST observations has changed dramatically since the advent\nof the satellite observations in 1982, i.e., the number of pre-1982 observations is remarkably small.\nIt is a well-known fact that the number of observations is fatally low, especially in the Southern\nOcean and in the deep winter season in high latitudes, and that there are not enough observations\nto fully resolve the SST seasonal cycle (see, for example, the spatial coverage and number of\nobservations archived in the World Ocean Database). I have doubts about the inclusion of the pre-\n1982 period in the analysis of this study based on some results read from the figures in the\nmanuscript, such as\n-Figure S1c-f: The winter trends in each hemisphere, which are important for defining the\namplitudes of the seasonal cycle, differ between products.\n-Figure S3: The magnitude and sign of the trend in the Southern Ocean varies greatly when the\ntime period is changed (1960- or 1982-) for both ERSST and HadISST (as also the authors mention\nin L103 but assume as minor). Note that the consistency between the two products does not resolve\nthis concern, as the two products use similar data sources prior to 1982, just with different mapping\nmethods.\n-Figure S1 a and b: There appears to be a large discrepancy between the model and observations\nin Figure S1 a and b (for a fair comparison, please use the same y-axis range for these panels).\nThis seems to be due to the strong questionable trend in the Indian Ocean sector of the Southern\nOcean that appears in the calculation since 1960. It is true that there could be a bias in the model,\nbut if the observations contain errors, wouldn't it be better to start from 1982 so that the model and\nobservations are more consistent?\n-Figure 1c: There seems to be no significant trend before 1980.\nResponse: Thank you for raising this issue. We acknowledge the limitations associated with pre-\nsatellite era SST observations, particularly in regions such as the Southern Ocean and during the\nwinter season at high latitudes. Taking your suggestion, we have decided to focus our analysis on\nsatellite-era data from 1982 to 2022 in the revised manuscript.\nWhile there are differences in the trend patterns between observations and simulations, we are\npleased to note that our main findings regarding the intensification of SST seasonality in the\nnorthern subpolar gyres and the northern flank of the ACC remain consistent across all datasets.\nWe thank you for highlighting this important issue, which helped to improve the robustness of our\nresults.\n(2) Conclusion about the mechanism\nThe authors concluded that the decrease in the annual mean MLD is the main \"mechanism\" for the\nincrease in the amplitude of the SST seasonal cycle by showing only an intermodel \"correlation\"\nbetween changes in the annual mean MLD and SST amplification. However, this is only for the\nglobal mean trend, and the spatial correspondence of the signals is not very good, although the\nauthors say it is good (comparing Fig. 3c to 1e). Even in the zonal mean (Fig. 3d), there is no\ndominant mixed layer thinning north of 50N, where the seasonal amplification is greatest, and the\nseasonal SST amplitude is not enhanced despite a strong mixed layer shallowing trend at 30S.\nAgain, this may be due to errors in the data product (in both IAP and Ishii) due to the lack of\nvertical profiles prior to the Argo era, as the model appears to be better (Fig. 3f).\nSince there may be several other mechanisms for the seasonal amplification of the SST besides\nthe one rejected by the authors (changes in the seasonal amplitude of the sea surface flux) (e.g.\nchanges in the annual mean sea surface heat flux and changes in the seasonal mean MLD), I believe\nthat the authors can only point out the relationship between the annual mean MLD and the SST\namplification, and not determine the mechanism (as the authors also mention in L203).\nResponse: We appreciate the insightful comments. We have made significant revisions to address\nthe concerns raised about the conclusion regarding the underlying mechanisms.\nFirst, we would like to address the difference between observed and simulated zonal mean MLD\nin the 50°-60°N region. This difference has been substantially reduced after presenting the trends\nin zonal-mean MLD over the period 1983-2022 (Fig. 5d, purple line). The previous analysis, based\non the average over 2013-2022, may have been strongly influenced by decadal variability.\nRecognizing the limitations of observational datasets, particularly in regions such as the Southern\nOcean, we have primarily relied on CMIP6 simulations in the revised manuscript, as these models\ncan capture the overall magnitude and spatial pattern of the intensification trend in key regions of\ninterest (Fig. 1), and the need to perform a detailed heat budget analysis.\nTo delve deeper into the underlying mechanisms, we conducted a careful mixed layer heat budget\nanalysis instead of the correlation analysis in the previous version. This approach allows us to\nquantitatively assess the contributions of various factors, including the thermal forcing term 𝑄\ninfluenced by the net surface heat flux and MLD, horizontal advection, and other residual\nprocesses (see Methods for more details). Through this approach, we identified decreasing MLD\nas the primary driver over most of the global ocean during the period 1982-2022.\nSpecifically, during 1982-2022, 𝑄 exhibits positive trends over most of the global ocean (Fig.\nR5b), particularly in the northern subpolar gyres. The residual term emerges as the major\ncontributor to the seasonal intensification in the equatorial oceans and also plays a role in the\nintensification on the northern flank of the ACC (Fig. R5d). In contrast, the contribution of\nhorizontal advection is considered negligible (Fig. R5c).\nFurther investigation is warranted into the primary factors shaping 𝑄 . These factors include\ndecreased MLD (Fig. R6a), reduced seasonal amplitude of MLD (Fig. R6b), increased heat uptake\nin the surface ocean (Fig. R6c), and enhanced amplitude of SHF over most of the global ocean\n(Fig. R6d). Through detailed decomposition, our analysis highlights the predominant role of the\nshallower MLD (Fig. R6e). Changes in the annual mean and seasonal amplitude of SHF also play\na role in the regional intensification of the SST seasonal cycle (Fig. R6c, d).\nIn the revision, we have rewritten the section “Mechanisms of the intensified SST seasonal cycle”\nto incorporate these important findings (Lines 175-222), which provide a compelling line of\nevidence for the dominant role of the decreasing MLD. Thank you for your valuable suggestions.\n(3) Possibilities for modeled dissolved oxygen\nIn general, even the Earth system models (ESMs) in the CMIP6 generation suffer from serious\ndissolved oxygen biases mainly due to incomplete parameterizations of air-sea exchange process\nand biological processes. Although it is difficult to obtain the observational trend of monthly sea\nsurface dissolved oxygen due to lack of historical observational data, the model climatological\nseasonal cycle should be validated using recent climatological observational data products, for\nexample, World Ocean Atlas (https://www.ncei.noaa.gov/products/world-ocean-atlas) and\nGridded Ocean Biogeochemistry from Artificial Intelligence – Oxygen\n(https://www.ncei.noaa.gov/access/metadata/landing-page/bin/iso?id=gov.noaa.nodc:0259304).\nResponse: Thank you for your comment. We have conducted a comparison of the climatological\namplitude of the dissolved oxygen seasonal cycle in the WOA18 and CMIP6. Our analysis,\ndemonstrates generally good agreement in terms of the spatial patterns and\nmagnitudes.\nIn the revised manuscript, we have updated Fig. 5a to include the seasonal cycle of dissolved\noxygen averaged over the Northern Hemisphere, represented by the dotted gray line.\nOther important concerns:\nL34: In references 8-11, I could not find any reference to “the effect of changes in the SST seasonal cycle\non marine ecosystems\". Please elaborate this sentence.\nResponse: The following literature has been cited:\nMueter, F. J., C. Broms, K. F. Drinkwater, K. D. Friedland, J. A. Hare, G. L. Hunt, W. Melle,\nand M. Taylor, 2009: Ecosystem responses to recent oceanographic variability in high-\nlatitude Northern Hemisphere ecosystems. Progress in Oceanography, 81, 93–110,\nhttps://doi.org/10.1016/j.pocean.2009.04.018.\nOttersen, G., S. Kim, G. Huse, J. J. Polovina, and N. C. Stenseth, 2010: Major pathways by\nwhich climate may force marine fish populations. Journal of Marine Systems, 79, 343–360,\nhttps://doi.org/10.1016/j.jmarsys.2008.12.013.\nL82：What does the plus/minus range mean? Please clarify the statistical methodology in the text.\nResponse: It represents the 5-95% confidence level of the linear trend, indicating the range within\nwhich 90% of the data falls. We have made it clear in the revised text (Line 82).\nL93: This is true only for the case of the Northern hemisphere (Fig. 1a and b). It seems that there\nare some differences in the case of the Southern hemisphere (Fig. S3a and b)\nResponse: Thank you for your comment. We have noted the discrepancy between the observed\nand simulated amplitude changes in the SH, particularly in the western South Pacific and the\nSouthern Ocean around 50oS. The lower pattern correlation in the SH (0.24 compared to 0.72 in\nthe NH) suggests the influence of large observational uncertainties, particularly in the data-sparse\nSouthern Ocean. Other factors, such as the low-frequency internal variability within the climate\nsystem and model deficiencies, may also have contributed to these differences. The difference\nbetween OBS and EXT is smaller than the value of the ensemble standard deviation over the\nmajority of global oceans.\nWe have highlighted this interhemispheric contrast in the revised text (Lines 120-124). Thank you\nfor pointing out this important issue.\nL113: How did the authors estimate the degree of freedom for the statistical test? Using the number\nof data gird point as the degree of freedom is overestimate because adjacent grids are not\nnecessarily independent. Please clarify the statistical method.\nResponse: Thank you. An “effective sample size” 𝑁∗ (Bretherton et al. 1999) is calculated and\nused to estimate the degrees of freedom in the significance test of correlation.\nWe have added a statement on the statistical significance test in the Methods section (Lines 383-\n387).\nL116 “Some discrepancies between observations and models may arise either from the influence\nof low-frequency internal variability”: Could you estimate the range of internal variability\n(ensemble standard deviation) using the ensemble members and compare the range with the\ndiscrepancies between observations and models?\nResponse: Thank you for your suggestion. In response, we have estimated the internal variability\nas the ensemble standard deviation of CMIP6 and other LENS ensemble members, respectively.\nAll the ensembles exhibit significant internal variability in regions where the seasonal\nintensification is prominent, particularly in the subpolar North Atlantic and Pacific.\nWe found that the difference between OBS and EXT is smaller than the value of the ensemble\nstandard deviation over most of the global ocean. This suggests that some of the discrepancies\nbetween observations and models may be internally generated by the climate system.\nWe have incorporated this discussion into the revised text (Lines 124-126).\nL215 “Another intriguing feature is the widespread amplitude reduction beneath the mixed layer,\nattributed to the overall shallowing of the mixed layer that impedes the propagation of the seasonal\ncycle signal into deeper ocean layers.”: Interesting.\nResponse: Thank you. We are currently conducting further investigations to understand this\nphenomenon in observations. We are excited about the potential insights this may provide and look\nforward to sharing our findings in future research.\nL242: “deoxygenation” rather than “oxygenation”?\nResponse: We apologize for the confusion, and have changed the term to \"oxygen solubility\" in\nthe revised text accordingly.\nL243 “The enhanced seasonal contrast of upper ocean oxygen level may foster the occurrence of\nhypoxic conditions characterized by exceedingly low oxygen concentrations by superimposing\nupon the long-term oxygen depletion trend over the past 60 years”: To quantify this impact, please\nshow the spatial map of the relative changes in the sea surface oxygen amplitude (i.e., oxygen\namplitude changes over the past 60 years divided by climatological seasonal amplitude).\nResponse: Thank you for your suggestion. The strongest intensification appears to be in the\nsubpolar North Atlantic (40-60oN, 60oW-0), with a relative change of over 15%. Intensification in\nthe subpolar North Pacific (40-60oN, 140oE-120oW) and Southern Ocean (35-50oS) also reaches\n5.5% and 6.9%, respectively. Note that the amplitude changes in the tropics (within 15o) are not\nshown because the small amplitude of the climatological annual cycle in these regions would generate\nextremely large ratio changes. Therefore, we have refrained from showing the ratio in the revision.\nHowever, we have incorporated the discussion into the revised text (Lines 259-262).\nL325：”This calculation is based on the climatology of the period 1970-2000 from the ERSST”:\nI could not get this. Please elaborate this.\nResponse: To clarify, we obtain the months of maximum and minimum SST based on the\nclimatological period from 1970 to 2000 using the ERSST dataset, and these months serve as the\nbasis for calculating the amplitude of the seasonal cycle for each year in both observations and\nsimulations.\nThis clarification has been included in the revised text (Lines 343-346)\nL338: “The MLD derived with the 0.01 kg m-3 threshold aligns perfectly with regions that exhibit\nan intensified seasonal cycle of upper ocean temperature at all latitudes”: I think this is result and\nthus this is not placed in the “Method”.\nResponse: Thank you, and we have reported this result in the text (Lines 225-226).\nFigure 1a and b: The left y-axis should be “Monthly trend” rather than “Amplitude trend”?\nResponse: Corrected as suggested.\nFigure S6 and Method: Why did not the authors use observational heat flux product\n(OAFlux https://climatedataguide.ucar.edu/climate-data/oaflux-objectively-analyzed-air-sea-\nfluxes-global-oceans, and JOFURO3 https://www.j-ofuro.com/en/)? Observational data is of\nlimited duration, but should be used for validation of atmospheric reanalysis.\nResponse: We appreciate your suggestion to use observational heat flux products for validation\npurposes. However, upon careful examination, we found significant uncertainties in the estimation\nof SHF products, including discrepancies in the sign of the changes. These\nuncertainties in SHF estimation have been reported in previous studies (e.g., Large & Yeager, 2009;\nYu et al., 2019), which makes it difficult to detect the long-term trends in the annual mean, let\nalone the amplitude of the SHF seasonal cycle.\nGiven the substantial uncertainties in the observational and reanalysis SHF datasets, and\nconsidering that climate models effectively capture the overall magnitude and spatial pattern of\nthe intensification trend in the SST seasonal cycle, we rely primarily on CMIP6 simulations and a\nheat budget-based attribution analysis to reveal the contribution of SHF change to the\nintensification of the SST seasonal cycle in the revision (Lines 217-222).\nFigure S6b: This analysis also should be done for each hemisphere like Fig. S7.\nResponse: Thank you for your suggestion. In response to your comment, we have conducted a\nsimilar analysis for each hemisphere, which shows a closer relationship between\nintensified amplitude in SHF and SST in the NH.\nHowever, as we addressed in your major comment #2, correlation does not imply causality.\nTherefore, we have employed an attribution method to quantify the specific contribution of SHF\namplitude changes to the intensified SST seasonal cycle. The results reveal that SHF amplitude\nchanges have contributions in the North Pacific and some modest contributions in the southern\nsubtropics. We have included this analysis in the revised manuscript (Lines 217-222).\nL176: What about the spatial correspondence between the heat flux changes and SST amplitude\nchange?\nResponse: Instead of conducting a correlation analysis, we have employed a mixed-layer heat\nbudget approach to directly evaluate the relationship between heat flux changes and SST amplitude\nchanges. While the primary contribution of SHF amplitude changes is evident in the North Pacific,\nit is relatively small compared to the dominant role played by the decreasing MLD.\nData: Is there any reason why the authors did not use the CESM2 large ensemble\n(https://esd.copernicus.org/articles/12/1393/2021/)? The CESM2 is one of CMIP6 generation\nmodels and it has outputs from single forcing experiment\n(https://journals.ametsoc.org/view/journals/clim/aop/JCLI-D-22-0666.1/JCLI-D-22-0666.1.xml).\nResponse: The reason we did not include the CESM2 LENS data is that their future projection\nsimulations are based on the SSP370 scenario rather than the SSP585 scenario, which we focused\non for consistency in our analysis. However, it is worth noting that CESM2 LENS also shows\nseasonal amplitude intensification with a similar spatial pattern to that in EXT.\nMethod: 0.01 kg m-3 is too small for the threshold when determining the MLD from the Grid data\nproduct. As noted in the reference, 0.01 kg m-3 is a value to use for raw CTD profiles, rather 0.1\nkg m-3 or 0.125 kg m-3 is more common because vertical profiles in gridded product are vertically\nsmoothed.\nResponse: Thank you for your comments. We have considered various criteria for calculating the\nMLD, including the commonly used thresholds of 0.125 kg m-3 and 0.03 kg m-3, as recommended\nby Montégut et al. (2004) for defining MLD from vertical density profiles.\nMLD estimates derived using the 0.03/0.125 kg m-3 criteria exhibit more pronounced increases in\nthe southern subtropics and the southern flank of the ACC. However, the decreasing\ntrends in MLD in the northern subpolar gyres and the northern flank of the ACC remain consistent\nunder all three criteria, highlighting the robustness of the shallower MLD estimates\nin these regions.\nWe have chosen the 0.01 kg m-3 criterion due to its alignment with the intensification of upper\nocean temperature seasonality, suggesting that most of the seasonal cycle\nchanges occur at this threshold. In contrast, the 0.125 kg m-3 criterion would yield significantly\ndeeper MLD estimates in the high-latitude oceans, where the seasonal cycle\nchanges are relatively small.\nFigure S4a: It would be better to show the spatial map of the CMIP6 model spread (standard\ndeviation) next to Fig. S4a.\nResponse: Thank you, and the standard deviation of the CMIP6 ensembles has been added as\nSupplementary Fig. 4b."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Thank you very much for addressing my quesfions and comments.\nI am very happy to read the new and more insighfful discussions.\nNow, this work can be accepted and published by Nature Communicafions.", "reviewer_2": "Second review of Human-induced intensified seasonal cycle of sea surface temperature by Liu et al.\nThank you for addressing my concerns and for your sincere responses to the comments. I believe that\nthe revised manuscript is much improved. In parficular, the inclusion of the heat budget analysis\nperformed by the authors has greatly clarified the discussion of the mechanism for the amplificafion of\nsea surface temperature. After addressing a few minor comments below, I can recommend this\nmanuscript to the journal for publicafion.\nMinor comments:\nL224 and L228: Fig. 5b, d?\nL281: The number seems not to be revised.\nFig. 1 capfion: OBS is the mean of ERSST, HadISST, and OISST?\nFig. 4: Please clearly menfion that this figure is a result from models (not observafion) in the capfion.\nL373: When the authors calculate amplitudes of fime-integrated tendency terms on the right-hand side\nof the equafion, they use values on the months of maximum and minimum of Tm (not of fime-integrated\ntendency terms), right? If so, please clarify it.\nSupplementary Fig. 1a: While the capfion says “OBS”, the panel shows “OISST”.\nThis review was also a valuable learning experience for me. Thank you.\nReviewed by Ryohei Yamaguchi"}, "author_response": "Thank you for your feedback and recommendation. We have addressed your comments to ensure\nclarity throughout the text.\nMinor comments:\nL224 and L228: Fig. 5b, d?\nResponse: Revised.\nL281: The number seems not to be revised.\nResponse: Revised.\nFig. 1 caption: OBS is the mean of ERSST, HadISST, and OISST?\nResponse: Revised.\nFig. 4: Please clearly mention that this figure is a result from models (not observation) in the\ncaption.\nResponse: Thank you, and we have revised as suggested.\nL373: When the authors calculate amplitudes of time-integrated tendency terms on the right-\nhand side of the equation, they use values on the months of maximum and minimum of Tm (not\nof time-integrated tendency terms), right? If so, please clarify it.\nResponse: Yes, thank you for pointing this out. We have made it clear in the revision (Lines 374-\n376).\nSupplementary Fig. 1a: While the caption says “OBS”, the panel shows “OISST”.\nResponse: The title of Supplementary Fig. 1a has been corrected."}]

Decision letter