Phase variations of the summer and winter seasons in the Bohai Sea during the last four decades

In most coastal oceans, the impacts of global warming on season duration and timing of seasonal transitions remain unknown. To mirror the reality of the ongoing climate change, the summer and winter seasons are redefined using the local water temperature thresholds in the Bohai Sea. Then the phase variations of these seasons are quantified using the duration and transition timing indices, including the duration (DUR), onset (ONS), and withdrawal (WIT) indices derived from the OSTIA SST dataset at a very high resolution (0.05°). During the last four decades (1982–2019), secular trends of summer indices extracted by the ensemble empirical mode decomposition (EEMD) method reveal that the summer DUR has an accumulated increase of about 17 days (4.5 days decade-1), which is primarily induced by the phase advance of the summer ONS by about 16 days (4.2 days decade-1). Spatial features of the duration and timing indices demonstrate that the lengthening of summer DUR and the phase advance of summer ONS have significantly enhanced in the shallow regions, due to the limited thermal inertia and the shorter period of the ocean’s memory. In contrast, the secular trend of winter DUR exhibits an accumulated shortening of about 18 days (4.8 days decade-1), which is induced by a moderately delayed winter ONS of 6 days (1.6 days decade-1) and a significantly advanced winter WIT of 12 days (3.2 days decade-1). The potential linkage between the phase variations in the oceanic seasonal cycle and those of the atmospheric forcing was investigated by analyzing both the interannual variability and the secular trend. Over the analysis period, the secular trend of an earlier summer ONS is related to a total reduction of cloud cover by 30% of its climatological mean and an increase of incoming solar radiation of 10 W m-2 month-1 in the late spring. Thus, our results highlight the influence of cloud cover in addition to wind speed on the temporal variations of season transition timing.

In most coastal oceans, the impacts of global warming on season duration and timing of seasonal transitions remain unknown. To mirror the reality of the ongoing climate change, the summer and winter seasons are redefined using the local water temperature thresholds in the Bohai Sea. Then the phase variations of these seasons are quantified using the duration and transition timing indices, including the duration (DUR), onset (ONS), and withdrawal (WIT) indices derived from the OSTIA SST dataset at a very high resolution (0.05°). During the last four decades , secular trends of summer indices extracted by the ensemble empirical mode decomposition (EEMD) method reveal that the summer DUR has an accumulated increase of about 17 days (4.5 days decade -1 ), which is primarily induced by the phase advance of the summer ONS by about 16 days (4.2 days decade -1 ). Spatial features of the duration and timing indices demonstrate that the lengthening of summer DUR and the phase advance of summer ONS have significantly enhanced in the shallow regions, due to the limited thermal inertia and the shorter period of the ocean's memory. In contrast, the secular trend of winter DUR exhibits an accumulated shortening of about 18 days (4.8 days decade -1 ), which is induced by a moderately delayed winter ONS of 6 days (1.6 days decade -1 ) and a significantly advanced winter WIT of 12 days (3.2 days decade -1 ). The potential linkage between the phase variations in the oceanic seasonal cycle and those of the atmospheric forcing was investigated by analyzing both the interannual variability and the secular trend. Over the analysis period, the secular trend of an earlier summer ONS is related to a total reduction of cloud cover by 30% of its climatological mean and an increase of incoming solar radiation of 10 W m -2 month -1 in the late spring. Thus, our results highlight the influence of cloud cover in addition to wind speed on the temporal variations of season transition timing. KEYWORDS phase variation, seasonal cycle, sea surface temperature (SST), secular trend, ensemble empirical mode decomposition (EEMD), cloud cover, Bohai Sea

Introduction
Sea surface temperature (SST) is an essential ocean variable for quantifying and understanding the ongoing climatic change. Seasonal cycle accounts for more than 90% of the total SST variance and thus has been regarded as the dominant variability of SST over midlatitudes (Dwyer et al., 2012). On a global scale, Thomson (1995) and Mann and Park (1996) proposed that anthropogenic climate change not only alters the amplitude of the seasonal cycle represented by the surface temperature, but also the phase of the season transitions. Although ecological consequences of the phenological mismatch in marine ecosystem, such as declined species recruitment and biodiversity loss, have been widely reported (Bertram et al., 2001;Asch et al., 2019;Staudinger et al., 2019;Visser and Gienapp, 2019;Laurel et al., 2021;Wilson et al., 2021), very few studies addressed the phase variations of the seasonal cycle of seawater temperature (Thomas et al., 2017;Alexander et al., 2018).
The Bohai Sea is a continental shelf sea with an average water depth of 18 m off the Northwest Pacific. With a total area of 77,000 km 2 , the Bohai Sea is semi-enclosed by one of the most urbanized areas of China involving several largest and densely-populated cities. Due to the limited heat capacity and environmental carrying capacity, the physical environment and ecosystem of the Bohai Sea are sensitive and vulnerable to the global and regional climate change as well as the highintensity human activities (Bian et al., 2016;Chen et al., 2019;Song and Duan, 2019;Wang et al., 2019). During 1970-2015, the significant warming of seawater temperature in the Bohai Sea caused an increased prevalence of harmful algal blooms at a rate of 40% decade -1 , with an earlier timing by 5.5 days decade -1 (Xiao et al., 2019). Since the seawater temperature determines the decomposition rate of organic matter, the onset and duration of the summer season also have potential influence on the frequency and intensity of seasonal hypoxia (Zhai et al., 2019). However, the temporal variations of the timing and duration of the seasons during recent decades are completely unclear at the regional scale.
Therefore, more investigation on the phase variations of the season transitions is urgently needed for the Bohai Sea. As an appropriate variable representing the season transitions, SST has been well measured and documented, especially over the satellite era (Merchant et al., 2019;Minnett et al., 2019). And when the water column is mixed seasonally, the changes in SST can be further transmitted to deeper waters (Alexander et al., 2018). Many recent studies have revealed that the SST time series and the associated length time series of the seasons are non-linear and non-stationary under dramatic climate change (Park et al., 2018;Xu et al., 2021;Wang et al., 2021a). In various climate change applications, the ensemble empirical mode decomposition (EEMD) method has been an objective and insightful means in revealing the secular trends from the non-stationary time series, including the timing and length time series of each season over drylands of the Northern Hemisphere midlatitudes (Wu et al., 2011;Franzke, 2012;Song et al., 2020;Wang et al., 2021a). Thus, we conducted the EEMD method to extract the secular trends of the lengths and timing time series of the summer and winter seasons in the Bohai Sea for further investigation on the phase variations of the seasonal cycle.
In this study, we provide the first quantitative assessment of the observed changes in both the lengths and timing of the summer and winter seasons in the Bohai Sea during the last four decades. Furthermore, the potential linkage between the phase variations in the oceanic seasonal cycle and those of the atmospheric forcing will be investigated by analyzing both the interannual variability and the secular trend.
2 Data and validation 2.1 Observational and reanalysis data 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.

Study area and validation
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  deviations. As a reference, Figure 2D shows the total number of insitu 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.

Definition of the seasons
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.

EEMD method
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).

Results
In this section, we first present climatology of the redefined summer and winter seasons, spatial structures of the summer/ winter ONS and WIT in the Bohai Sea, and interannual variability and secular trends of the summer/winter DUR, ONS and WIT during 1982-2019. Subsequently, to investigate the possible role of physical drivers on the observed phase variations of the summer and winter seasons, we examined the potential relationship between the associated length and timing indices and the atmospheric variables, including the net heat flux, all the heat flux components (solar and longwave radiation, latent and sensible heat flux), zonal and meridional wind components, wind speed and cloud cover.

Redivision of the summer and winter seasons
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  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.

Secular trends of length and timing
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.

Summer ONS and cloud cover
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 (spatiallyaveraged 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.

Winter ONS and wind speed
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.

Winter/summer WIT and surface heat flux
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  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.

Summer expansion and marine heatwaves
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.

Local phase variation and global warming
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 Correlation coefficients between time series of (A) the late fall turbulent heat fluxes and the winter ONS, (B) the late fall wind speed and the winter ONS. Only significant correlations at the confidence level of 95% are shown in the ice-free regions of the Bohai Sea. Black contour is the same as 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.

Summary and conclusions
Understanding the phase variations of the summer and winter seasons is critical to achieving effective marine conservation and making appropriate adaption policies for the coastal ocean regions. In the Bohai Sea, both the duration and timing of the climatic summer and winter seasons have changed substantially during the last four decades . Using the EEMD method, secular trends of summer/winter DUR reveal that the summer season has lengthened by17 days (4.5 days decade -1 ) over the entire Bohai Sea, whereas the winter season has shortened by 18 days (4.8 days decade -1 ) in the ice-free regions.
In the Bohai Sea, the lengthening of summer DUR is primarily induced by the significant advance of summer ONS (4.2 days decade -1 ), rather than the very slight delay of summer WIT. Both the secular trend and inter-annual variations of the summer ONS are highly influenced by those of the late spring solar radiation and heat gain at the sea surface, which are significantly determined by the amount of cloud cover over the same period. In the ice-free regions, the temporal variations of winter ONS/WIT are significantly influenced by those of the late fall/winter turbulent heat fluxes at the sea surface, which are related to the magnitudes of local wind speed over the same period. Therefore, the atmosphere circulation can modulate the duration and transition timing of the climatic summer and winter seasons in the coastal oceans through various thermal and dynamic processes.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.