Daily mean SST during 1982–2019 was obtained from the operational sea surface temperature and sea ice analysis (OSTIA, Donlon et al., 2012; Good et al., 2020) dataset developed by the U.K. Meteorological Office and distributed through the Copernicus Marine Environment Monitoring Service (CMEMS). The OSTIA dataset has a very high spatial resolution of 0.05° thus provides the advantage to resolve the coastal ocean processes in the marginal seas (Minnett et al., 2019). In this study, the OSTIA SST data has been assessed using in-situ observations measured at fixed coastal stations operated and maintained by the State Oceanic Administration of China.

Daily mean fields of wind velocity, air temperature, total cloud cover (TCC), surface heat flux components including solar radiation, longwave radiation, latent and sensible heat fluxes during the same analysis period were obtained from the European Centre for Medium-Range Weather Forecast (ECMWF) reanalysis dataset ERA5 at 0.25° of horizontal resolution. The dataset for marine heatwaves (MHWs) intensities was obtained from the database Science Data Bank (Zhang et al., 2022; http://www.scidb.cn/en/s/nqauYn). The daily intensities of MHWs are calculated based on the difference between the National Oceanic and Atmospheric Administration Optimum Interpolation SST and its 30-year climatological mean at 0.25° grid.

The Bohai Sea ( Figure 1 ) can be divided into three coastal bays named Laizhou Bay, Bohai Bay and Liaodong Bay, a shallow basin named the Central Bohai Sea and the Bohai Strait. At Laizhou Bay, the most area of Bohai Bay and the head of Liaodong Bay, water depth is less than 18 m, which is the averaged water depth of the entire Bohai Sea. Then the water depth gradually deepens from these coastal bays to the Central Bohai Sea and the Bohai Strait. The in-situ observation stations of Beihuangcheng (BHC), Qinghuangdao (QHD) and Bayuquan (BYQ) (shown in Figure 1 ) provide the longest running oceanographic time series in the study areas.

The Taylor diagram (Taylor, 2001) can provide a graphical representation of how close the OSTIA SST matches the in-situ SST observation in the Bohai Sea during the satellite era. Figures 2A–C show that the correlations between the OSTIA and the in-situ observed SSTs are between 0.95 and 0.99 at the stations of BHC, QHD and BYQ, which indicate the OSTIA SST time series can closely reflect the temporal variability derived from the in-situ observations. The in-situ observed SST has been stored as the 10-day average record. Due to the differences in sampling intervals, the remote sensing SST has more fluctuations which caused the RMS deviations. As a reference, Figure 2D shows the total number of in-situ observation records in each year for all the oceanographic stations. The most consecutive missing records occurred during 1985–1990 at the stations of BHC and QHD.

Traditionally, both the meteorological and oceanographic seasons divide the year into four three-month periods. In the meteorological seasons, as being the warmest quarter of the year, summer is defined as June, July, and August; And as being the coldest quarter, winter is defined as December, January, and February at midlatitudes (Trenberth, 1983). The oceanographic seasons correspondingly postpone the meteorological seasons by one month due to its higher heat capacity and larger thermal inertia. However, the traditional definition of oceanographic seasons with fixed length and onset date does not fully agree with the recent observed climate extremes (e.g., marine heatwaves) and the detected phase variations in the seasonal hydrological cycle (Klein Tank and Können 2003; Thomas et al., 2017). Thus, the oceanographic seasons in the coastal oceans are redefined using the local SST to provide a more appropriate and objective perspective under climate change.

In this article, we use the local seawater temperature criteria to redefine the summer and winter seasons in the Bohai Sea. The local temperature thresholds for the summer and winter seasons are defined as the 75th and 25th percentile temperature values obtained from the climatological annual cycle of the daily water temperature averaged over the period of 1982–2019, which indicate the warmest and coldest quarters in each year, respectively. Applying the method of Christidis et al. (2007) and Park et al. (2018), summer onset (ONS) and withdrawal (WIT) are defined as the calendar dates when the local SST starts to exceed and decreases to below the corresponding threshold, respectively; Winter timing indices are defined following the same logic. And summer/winter duration (DUR) is defined as the number of calendar days from ONS to WIT. To avoid any ambiguity, the winter of 1982 refers to 1982/83 winter in this study.

As an important parameter in this season redivision algorithm, the temperature threshold is calculated on the basis of local temperature throughout the whole year during 1982–2019. It is noted that the winter temperature under the sea ice has been relaxed to –1.8°C at a complicated timescale related to the ice concentration as described in section 2.4 in Good et al. (2020). Thus, we perform the sensitivity tests on the threshold values to verify the rationality of redivision method and threshold setting. First, we compare the summer season indices including ONS, WIT and DUR calculated using the original threshold of 75th percentile with those using the thresholds of the original value minus/plus 0.5°C and 1.0°C, respectively. These parallel tests clearly reveal that both the interannual variations and secular trends of the summer season indices calculated with the slightly different thresholds are substantially similar in the Bohai Sea, although their actual values are somewhat different (see Figure S1 in the supplemental material). Therefore, the phase variations of the summer season indices are insensitive to the threshold changes in the magnitude that can be induced by sea ice. Further, we compare the summer season indices calculated using the original threshold of 75th percentile with those using different percentile thresholds of the 70th and 80th percentiles; Parallel tests for the winter season are performed using the thresholds of 20th, 25th and 30th percentiles, respectively. The temporal variability of both the summer and winter season indices does not show significant difference in the sensitivity tests with different thresholds (see Figures S2 and Figure S3 in the supplemental material). Thus, the applied threshold criteria are reasonable in this study and consistent with the idea that summer/winter is the warmest/coldest quarter of the year.

The ensemble empirical mode decomposition (EEMD) method can decompose the original non-stationary time series into a set of oscillatory components at different time scales and a secular trend better reflecting the underlying physics (Wu et al., 2011; Ji et al., 2014; Wang et al., 2021a). Over the Bohai Sea, time series of DUR, ONS and WIT of the summer and winter seasons are non-stationary of which statistical properties changing through time ( Figures 3B–G ). Therefore, we applied EEMD method to decompose the length and timing indices time series of the above two seasons to extract the oscillatory components as well as the secular trends over the last four decades. The oscillatory components, also known as the intrinsic mode functions, whose amplitude and frequency can vary over time, thus provided a more accurate time-frequency representation than the sine and cosine components in the traditional Fourier Transform (Stallone et al., 2020). And the EEMD secular trends were used in further attribution analysis of the phase variations in the seasonal cycle of water temperature.

From the spatial perspective, EEMD is a completely local method, for the decomposition of time series at one grid is independent of the data characteristics at other grids (Ji et al., 2014). Thus, spatial structures of the secular trends of ONS/WIT timing indices extracted by EEMD were obtained to investigate how the local water depth influences the phase of water temperature cycle in the Bohai Sea. In EEMD calculation, the added noise series has an amplitude of 0.2 times the standard deviation of the corresponding data and the ensemble number is set to 100 (Ji et al., 2014; Gaci, 2016).

In this study, the summer and winter seasons during 1982–2019 are redefined using local water temperature thresholds in the Bohai Sea. As summer criteria, 75th percentiles of SST ranged from 18.9°C in the Bohai Strait to 23.4°C at the top area of Laizhou Bay; Corresponding winter criteria of 25th percentiles of SST ranged from 1.2°C at the top area of Liaodong Bay to 5.5°C in the Bohai Strait. The localized temperature thresholds can provide more reasonable criteria for the season division. Based on these criteria, time series of the length and timing indices, including DUR, ONS and WIT were calculated for the summer and winter seasons, respectively.

The climatological length and timing indices are averaged over the entire Bohai Sea and illustrated in Table 1 and Figure 3A . On average, the summer season started on 25 June and ended on 27 September, while the winter season started on 27 December and ended on 29 March. As shown in Figure 4 , spatial distributions of climatological ONS and WIT for both the summer and winter seasons showed clear dependency on the water depths. Top areas of Laizhou Bay, Bohai Bay and Liaodong Bay, where water depths are less than 18 m, exhibit earlier ONS and WIT of about 10 days compared with those averaged over the Bohai Sea for the summer and winter seasons; On the contrary, the Bohai Strait with the deepest water depth shows the latest seasonal transitions with more than 10 days compared with the Bohai averages. These spatial features of timing indices indicate that the significant difference in thermal inertia determined by the local water depth can highly influence the timing of seasonal transition in the Bohai Sea.

Over the analysis period, time series of the length and timing indices for the summer and winter seasons exhibits significant inter-annual fluctuations ( Figures 3B–G ). For the summer season, the longest and shortest duration of 108 days and 78 days occurred in 2019 and 1985, respectively; For the winter season, the longest and shortest duration of 112 days and 72 days occurred in 2009 and 2016, respectively. During the last four decades, DUR time series for the summer and winter seasons exhibits significant lengthening and shortening trends, respectively. Since the traditional statistical method is not appropriate to extract the secular trends from these non-stationary time series, EEMD method will be applied to quantify the secular trends more accurately in the next session.

First, we conducted the EEMD method on the time series of the summer/winter length and timing indices and investigated the temporal characteristics of their secular trends during the period of last four decades. Figure 5 shows EEMD decompositions of the summer indices averaged over the entire Bohai Sea and those of the winter indices averaged over the ice-free areas as illustrated in the right panel of Figure 6 . For the summer season, the secular trend of DUR exhibits a continuous lengthening during 1982–2019, with an accumulated length increase of about 17 days over 38 years (4.5 days decade^-1). In line with the global warming, the summer ONS and WIT has shown monotonic advance and delay in their secular trends, respectively. Compared with the slightly delayed summer WIT, the summer ONS advanced by about 16 days (4.2 days decade^-1) contributes significantly more to the lengthening of the summer season during the last four decades. In contrast, the secular trend of winter DUR exhibits an accumulated shortening of about 18 days during 1982–2018 (4.8 days decade^-1), which is induced by a moderately delayed winter ONS of 6 days as well as a significantly advanced winter WIT of 12 days (lower panel in Figure 5 ). Hence, the secular trends derived from the EEMD analysis clearly demonstrate an extension of summer, a contraction of winter, and thus the significant phase variations in the seasonal cycle of SST in the Bohai Sea.

Second, we further investigated the spatial features of accumulated variations of the summer/winter length and timing indices over the last four decades at each grid point in the Bohai Sea by means of the spatial locality of EEMD method. Figure 6 illustrates spatial structures of the accumulated variations of DUR, ONS and WIT are derived from the corresponding secular trends grid by grid. For summer, the spatial patterns suggest that the long-term variation of summer DUR is basically consistent with that of the summer ONS. The top areas of Liaodong, Bohai and Laizhou bays exhibit the strongest lengthening trends (by 30 days) of the summer DUR and also the most significant advance (exceeding 20 days) of the summer ONS, whereas the Bohai Strait area shows no significant variation in the summer length and timing indices, due to the lateral heat transport with the northern Yellow Sea (left panel in Figure 6 ). Previous studies found that the thermal inertia of a water column on the continental shelf is linearly proportional to its local water depth (Huang et al., 1999; Xie et al., 2002; Bian et al., 2016). Thus, both the lengthening of summer DUR and the advance of summer ONS tends to be more significant in the shallow (< 18 m) regions, due to the smaller magnitude of the ocean’s heat capacity and the shorter period of ocean’s memory.

Since the SST in OSTIA dataset has been deliberately relaxed towards –1.8° C in the sea ice covered regions (Good et al., 2020), our analysis for the winter season mainly focuses on the ice-free regions of the Bohai Sea. During 1982–2018, the shortened winter DUR (by 20–30 days), delayed winter ONS (by 0–10 days) and earlier winter WIT (by 10–20 days) can be found in the deep and ice-free regions of the Bohai Sea, with the water depths spanning 18–86 m (right panel in Figure 6 ). In general, the secular trends of the season durations and their transition timing have the same positive/negative signs in the Bohai Sea, which indicates their common connections with the atmospheric forcing.

By definition, the timing of summer ONS is closely related to the water temperature rising process in late spring, which is primarily influenced by the heat gain due to the incoming solar radiation. During the last four decades, time series of the late spring solar radiation is obtained by accumulating the incoming solar radiation within a month prior to the summer ONS in each year. Figure 7A illustrates that time series of the late spring solar radiation is significantly correlated with that of the summer ONS over most of the Bohai Sea. The negative correlation (spatially-averaged r = –0.65) indicates that the larger late spring solar radiation leads to an earlier summer ONS. Furthermore, the cloud cover is regarded as an important regulator of the amount of incoming solar radiation that reaches the sea surface (Hughes, 1984; He et al., 2019; Gao et al., 2020; Kejna et al., 2021); And the cooling action of local cloud predominates over the entire Bohai and Yellow Seas (Yuan, 2011; Yao et al., 2020). Figure 7B shows that the significant positive correlations (r ranging from 0.5 to 0.8) between time series of the one-month accumulated TCC in late spring and that of the summer ONS exhibits in large areas of the Bohai Sea. Meanwhile, the spatial pattern of correlation coefficients between the late spring TCC and summer ONS is very similar to that of the coefficients between the late spring solar radiation and summer ONS ( Figures 7A, B ). Hence, it suggests that the reduced TCC amount in late spring could induce an increased solar radiation at the sea surface, thus enhancing the SST rising process and causing an earlier summer ONS.

We further explored the statistical relationship between the temporal evolution of the late spring cloud cover and that of the summer ONS over the entire Bohai Sea. During 1982–2019, time series of the late spring TCC has a negative correlation r = –0.87 with that of the incoming solar radiation at the sea surface and a positive correlation r = 0.73 with that of the summer ONS, respectively (right panel in Figure 7 ). From the perspective of long-term variations, the amount of late spring TCC over the Bohai Sea clearly shows a significant decreasing trend of 0.16 over the recent 38 years, which approximately accounts for 30% of the climatological mean TCC of 0.53 ( Figure 7C ). This decreasing trend of TCC is also confirmed by other reanalysis and ground station datasets (Zong et al., 2013). The dramatic loss of late spring TCC caused a long-term increase of 10 W m^-2 month^-1 in the magnitude of incoming solar radiation during the same period and further induced a phase advance of 16 days in the summer ONS ( Figures 7D, E ). In conclusion, the timing of summer ONS is most likely determined by the late spring cloud cover during the last four decades.

Changes in surface wind speed can potentially influence the SST decreasing process in late fall and the corresponded winter ONS through the sea surface heat fluxes. The significant positive correlations (r ranging from 0.6 to 0.8) between the time series of the winter ONS and that of the turbulent heat fluxes accumulated within two months prior to the winter ONS are found over the entire ice-free regions of the Bohai Sea ( Figure 8A ). These significant correlations suggest that the temporal variations of winter ONS are possibly controlled by that of the surface cooling process influenced by the turbulent heat fluxes. The negative correlations (r ranging from –0.8 to –0.5) further indicated that lower wind speed in late fall favors the weakened surface cooling, which results in a delayed winter ONS ( Figure 8B ).

We then investigated the statistical relationship between the temporal evolution of the wind speed and that of the winter ONS over the ice-free region of the Bohai Sea. As illustrated in Figures 8C–E , correlation analysis revealed that the time series of late fall wind speed is negatively correlated with both of the turbulent heat fluxes (r = –0.68) and winter ONS (r = –0.67), respectively. And the secular trends suggested that the weakened late fall wind speed has caused a decrease of 15 W m^-2 month^-1 in the magnitude of turbulent heat fluxes during the same period. In summary, both the correlation analysis and the secular trends derived from the EEMD analysis support that the weakened local wind speed in late fall could result in the decreased surface heat loss and the gradually postponed winter ONS.

Timing of the winter/summer WIT is determined by temperature rising/decreasing process over the period of the late winter/summer, which are influenced by the heat gain/loss at sea surface. By analyzing the time series of total heat flux and its associated components including solar radiation, longwave radiation and turbulent heat fluxes (including latent and sensible heat fluxes), we find that the temporal variation of turbulent heat fluxes is more responsible for that of the surface heat gain/loss. Thus, correlations between the time series of turbulent heat flux and the winter/summer WIT are investigated over the Bohai Sea. In Figure 9A , time series of the turbulent heat fluxes in late winter shows a moderate correlation (r = –0.66) with that of the winter WIT. Negative correlation indicates that the greater magnitude turbulent heat fluxes favors to a slower temperature rising process and a delayed winter WIT. In Figure 9B , the turbulent heat fluxes were further related to the wind speed over the same period, thus time series of the late winter wind speed shows a moderate correlation (r = 0.45) with that of the winter WIT. Thus, larger wind speed in the late winter favors to a delayed winter WIT. Similarly, temporal variation of the summer WIT is determined by the surface heat loss, which is closely related to the turbulent heat fluxes and wind speed in late summer. As illustrated in Figures 9C, D , the summer WIT is moderately correlated with both the turbulent heat fluxes (r = 0.57) and the wind speed (r = –0.55) accumulated over the late summer. Thus, larger wind speed in the late summer favors to a faster temperature decreasing and an earlier summer WIT.

From the perspective of the secular trend, the long-term variation of summer WIT during the last four decades is very limited ( Figure 5C ); Whereas the long-term variation of winter WIT shows a significant advance of 12 days ( Figure 5F ), induced by the total decrease of 21 W m^-2 month^-1 in the magnitude of the late winter turbulent heat fluxes during the analysis period of this research.

In the traditional meteorological/oceanographic seasons, there are four seasons, which encompass three months each and generally match the past climatological temperature patterns at midlatitudes. During the past few decades, the significant increase of the frequency and severity of ocean extreme events including MHW has been revealed and confirmed by both the long-term observations and the climate model simulations (Sen Gupta et al., 2020; Li and Donner, 2022). However, the fixed-length summer season does not show any linkage to the occurrence and intensity of MHW. Thus, a more reasonable definition of the seasons is needed to mirror the reality of recent climate change more closely.

In this study, we provide a temperature-based redivision of the summer and winter seasons. To further explore the potential relationship between this climatic summer season and the MHWs, time series of the summer DUR and accumulated MHW intensity during 1982–2019 are illustrated in Figure 10 . Here the accumulated MHW intensity is defined as the accumulated sum of daily intensity during all the MHW periods in each year, which represents the overall SST anomalies and heat stress from the extreme warm events. In the Bohai Sea, the accumulated MHW intensity has a moderate correlation (r = 0.58) with the summer DUR; Meanwhile, secular trend of the accumulated MHW intensity has a very high correlation (r = 0.94) with that of the summer DUR. It is indicated that the MHW intensity and the associated heat stress were significantly stronger in a lengthened climatic summer season.

Therefore, the temperature-based definition of the seasons can better reflect the climatic extreme events. And the research on the increased and enhanced MHWs should be conducted with full consideration of the ongoing and future oceanic summer expansion, which is also suggested by Wang et al. (2021b) and Park et al. (2022) in the heatwave study for the land areas.

During the past few decades, global warming leads to the significant summer expansion and winter contraction across the Northern Hemisphere midlatitudes; However, the magnitudes of length variations are not uniform on the spatial scale. Park et al. (2018) suggested that the summer season had lengthened from 12 to 36 days during 1953–2012 (2.0–6.0 days decades^-1) over the continental Eurasia; Wang et al. (2021a) reported that the lengthening/shortening rate (5.1/2.3 days decades^-1) of the summer/winter season was particularly enhanced in drylands, where the regional warming was greater than the global average. In this study, the secular trends of summer and winter lengths in the Bohai Sea were generally consistent with those hemispheric-scale trends. It is revealed that the regional lengthening rate of summer (4.5 days decades^-1) was close to that in the land area experiencing rapid warming, whereas the shortening rate of winter (4.8 days decades^-1) was approximately twice the land value.

Previous temperature and heat content budgets for the Bohai Sea indicate that their spatiotemporal variations are primarily controlled by the local air-sea heat flux and the associated atmospheric forcing, which can be further linked to large-scale climate variability. Hence, local phase variations of season transitions derived from the water temperature cycle are correspondingly attributed to the atmospheric factors including cloud cover and wind speed. During late spring, cloud formation and its vertical structures are highly related to the onset and northward advance of the East Asian Summer Monsoon, which is modulated by the Pacific and Indian Ocean SST, the snow cover on the Tibetan Plateau and the location of Western Pacific Subtropical High (Sun et al., 2019; Gao et al., 2020). During late fall, wind speed is primarily related to the East Asian Winter Monsoon, which is influenced by the East Asian Trough, Aleutian Low, North Pacific Oscillation and ENSO. In further quantitative research, global climate models can be important tools for the investigation on the linkage between the phase variations of the seasons in the coastal ocean and the large-scale climate variability under global warming.