In this study, the Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA) SST products with a 0.05° resolution, covering the 1985-2022 period in the Eastern China Seas, were downloaded from the Europe Copernicus program (http://marine.copernicus.eu). The cloud-free, global SST product provided by the global ocean data assimilation experiment High-Resolution SST Pilot Project (GHRSST-PP) merges OSTIA SST data from multi-sensor satellite remote sensing observations (Donlon et al., 2012). This has been previously validated against in situ observations within the study area, showing strong agreement with a correlation coefficient of 0.94 and a slight bias of 0.12°C, confirming its reliability and accuracy in local conditions and outperforming other SST products (Wang et al., 2020a).

Daily National Centers for Environmental Prediction (NCEP) meteorological products for 1985-2022 were obtained from the Asia-Pacific Data-Research Center (http://apdrc.soest.hawaii.edu/data/data.php). The collected parameters are air temperature at 2 m and latent heat flux (Q_LH). NCEP is a third-generation global coupling seasonal forecast reanalysis data with a spatial resolution of 0.3125° (~38 km). It adopted advanced global data assimilation system and database, carried out quality control and assimilation processing of observation data from various sources (ship, aircraft, satellite, etc.), and obtained a set of reanalysis data (He and Zhao, 2018; Zhu et al., 2021). NCEP’s daily meteorological products are widely used by meteorologists, climate scientists, and researchers for studying weather patterns, climate variability, and long-term climate trends (He and Zhao, 2018). Their accuracy and reliability make them invaluable tools for understanding and predicting weather phenomena on regional and global scales (He and Zhao, 2018; Zhu et al., 2021).

The numerical values of the Multivariate ENSO (El Niño/Southern Oscillation) Index Version 2 (MEI) were obtained from the Physical Sciences Division (PSD) of the Earth System Research Laboratory, part of the National Oceanic and Atmospheric Administration (NOAA) (https://psl.noaa.gov/enso/mei/). Negative MEI values indicate the cold ENSO phase (La Niña), whereas positive MEI values indicate the warm ENSO phase (El Niño; Wolter and Timlin, 1998).

In our analysis, several assumptions were made to ensure the robustness and validity of the results. We assumed that the data from different sources, such as the OSTIA and NCEP datasets, are homogeneous and compatible for combined analysis. This assumption is based on the standardization procedures and calibration techniques employed by these data providers, ensuring that their datasets are consistent over the study period. For missing data points, we employed a linear interpolation method, assuming that the changes between observed values are linear. This approach is commonly used in climatological studies and helps in maintaining the continuity of the time series (Sattari et al., 2017). We assumed that the underlying statistical properties of the climate variables, such as the mean and variance, remain constant over the period of study. This assumption is crucial for the validity of trend analysis and is typically justified in climatological studies over decadal timescales. These assumptions were carefully considered and tested to minimize potential biases in our analysis.

The SST in the Eastern China Seas exhibits seasonal characteristics, with the highest SST typically occurring in summer (Si et al., 2016). The specific day when SST reaches its peak in summer is varies. We extracted the SST time series for each pixel annually and smoothed with a 31-day moving average. Following the method illustrated in Figure 2 (where the maximum value in the time series is designated as T_max, and its corresponding time is T_max timing), we identified the T_max and T_max timing for each pixel each year. We examined the potential autocorrelation in T_max and T_max timing across different years and found no evidence of autocorrelation among a sample of pixels. This suggests that the T_max and T_max timing in any given year is independent of the T_max and T_max timing in preceding years. The decadal trends of T_max and T_max timing were estimated by fitting a linear regression relating T_max and T_max timing to year for each of the pixel’s yearly time series. The slopes of the fitted linear regressions provide an estimate of the rates of change for each pixel (units: °C decade^-1 for T_max, and days decade^-1 for T_max timing). To assess whether there were significant trends in the T_max and T_max timing time series, we calculated the Pearson correlation coefficient (r) between the original values and the values obtained from linear regression. A t-test coefficient p can be calculated to determine the significance level. A significant trend in the time series is considered to exist when p<0.05. We also utilized the r to measure the linear relationship between two continuous variables, which ranges from -1 to 1. A correlation coefficient of 1 indicates a perfect positive correlation, -1 indicates a perfect negative correlation, and 0 indicates no linear relationship. Similarly, significance levels are determined using a t-test, where p<0.05 indicates a significant correlation.

Similar to SST, air temperature also peak during the summer. We employed the same method as for SST to obtain the maximum values and corresponding times for air temperature (AirT_max and AirT_max timing).

This study follows the methodology described by Hobday et al. (2016) and Oliver et al. (2018). Marine heatwaves (MHWs) are identified at each pixel using daily OSTIA data. An MHW event is defined when the daily SST exceeds the 90th percentile threshold, calculated from a 38-year dataset of daily SSTs. To determine seasonally varying 90th percentile thresholds, daily SSTs are computed within an 11-day window centered on the corresponding date of each year from 1985 to 2022. Subsequently, a 31-day moving average is applied to smooth the 90th percentile thresholds. Any contiguous days with temperatures above the threshold for at least five consecutive days and a break of no more than three days is considered as a single event. Matlab codes we used in this study to calculate MHW characteristics are available at https://github.com/ZijieZhaoMMHW. The Proportion of MHWs at T_max timing is defined as the number of pixels in the study area where MHWs occur at T_max timing divided by the total number of pixels.

The increase in ocean temperature has gained significant interest due to its connection with global warming (Trenberth et al., 2007; Rhein et al., 2013). Numerous studies have reported a warming trend in summer SST in the Eastern China Seas, with a consensus that the rate of warming is significantly higher than the global average (Cai et al., 2017; Tang et al., 2020; Wang et al., 2023). T_max, representing the annual summer extreme temperature in the study areas, has been warming at a rate of 0.3°C decade^-1, which is 52% higher than the global average and 94% higher than the reported summer average warming rate in other studies (Cai et al., 2016, 2017; Wang et al., 2023). The notable increase in extreme temperatures poses a greater threat to the ecosystem, as organisms may already be operating at the upper limit of their thermal tolerance (Wernberg et al., 2013; Chaidez et al., 2017). Although the ability of organisms to adapt to warming through a poleward shift in their biogeographical range has been shown (Wernberg et al., 2016), the significant rise in T_max may have major impacts on the migration of fishery resources in the Eastern China Seas.

It can be observed from Figure 4 that there are significant differences between the southern and northern parts of the Eastern China Seas in terms of the T_max warming rate, change in T_max timing, and the relationship between T_max and T_max timing. We divided the study area into southern and northern parts using 31.5°N as the boundary and analyzed the changes in T_max and T_max timing in both the southern and northern parts separately, as shown in Figure 5 . In both the southern and northern parts, T_max increased significantly (p<0.01), with a warming rate of 0.4°C decade^-1 in the northern part, which is twice the rate of 0.2°C decade^-1 in the southern part. In the northern part, T_max timing significantly decreased over time (p=0.04), meaning the occurrence time of T_max advanced significantly, with a change rate of 1.3 days decade^-1, and T_max and T_max timing showed a significant negative correlation (r=-0.5, p<0.01). In the southern part, T_max timing significantly increased over time (p=0.04), meaning the occurrence time of T_max significantly delayed, with a change rate of 2.8 days decade^-1. The T_max in the northern region are lower than in the south ( Figure 3A ), and the faster increase in T_max in the north may lead to a reduction in the T_max difference between the north and south. T_max occurs earlier in the south than in the north ( Figure 3B ), and the significant increase in T_max timing in the south and the significant decrease in T_max timing in the north may lead to a reduction in the difference in T_max occurrence times between the north and south. Variance analysis shows that T_max (r=-0.54, p<0.01) and T_max timing (r=-0.36, p=0.03) in the study area are both significantly negatively correlated with T_max (Figure not show). Thus, while T_max is significantly warming, the difference in T_max and T_max timing between the southern and northern parts is diminishing.

The relationship between air temperature and SST is complex and interconnected (Al-Shehhi, 2022; He et al., 2023). Changes in air temperature can influence SST through various mechanisms, resulting in a strong correlation between air temperature and SST (Feng et al., 2018). In the Eastern China Seas, air temperature reaches its peak during the summer, approximately 6.4 days earlier than the peak SST in the southern part and 9.5 days earlier in the northern part. There is a high correlation between AirT_max and T_max, with the correlation coefficients of 0.93 in both the southern and northern parts, resulting in nearly identical trends for AirT_max and T_max in both parts ( Figure 6 ).

However, the correlation between AirT_max timing and T_max timing differs between the south and the north. In the southern part, AirT_max timing shows a similar increasing trend to T_max timing, whereas in the northern part, the difference between AirT_max timing and T_max timing has significantly decreased (p<0.01) ( Figure 6 ). This reduction in the time lag between T_max and AirT_max, along with the lack of a clear trend in AirT_max timing, results in an advancing trend in T_max timing in the northern part ( Figure 6 ).

Higher SST often leads to increased Q_LH due to enhanced evaporation rates (Kumar et al., 2017; Mohan et al., 2023). In the northern part of the Eastern China Seas, there is a strong correlation between Q_LH and T_max in August ( Table 1 ). This increased Q_LH can influence SST by promoting heat transfer from the ocean surface to the atmosphere, potentially moderating further increases in SST (Sekizawa et al., 2018; Ganguly et al., 2024). In August, strong Q_LH suppresses the continuous increase in SST in the northern part, causing T_max timing to advance. This may be an important reason for the significant negative correlation between Q_LH and T_max timing ( Table 1 ). In contrast, this relationship between Q_LH, T_max and T_max timing is much weaker in the southern part. In the southern part, approximately 39% of T_max occurs in July, and 61% occurs in August. The Q_LH in both July and August shows a significant increase in the south ( Figure 7 ). High Q_LH in July suppresses the increase in SST, delaying T_max to August, which diminishes the effect of high Q_LH in August on advancing T_max. In the northern part, 93% of T_max occurs in August, and there is no significant increase in Q_LH in July. This makes the impact of the increased Q_LH in August on advancing T_max particularly evident. Additionally, regional variations in Q_LH can significantly impact local SST trends, emphasizing the importance of understanding these dynamics for accurate climate modeling and forecasting in the Eastern China Seas.

The MEI is an important indicator of ENSO variability, which has been shown to influence SST patterns globally (Yuan et al., 2018; Cai et al., 2019). In the Eastern China Seas, ENSO events can significantly impact summer SST dynamics (Sun et al., 2021; Liu et al., 2023). There is a significant negative correlation between MEI and T_max during July and August in the southern part (Table 1), indicating that La Niña tends to favor higher T_max values. During El Niño and La Niña events, there are notable shifts in both the position and intensity of the western Pacific subtropical high (WPSH) (Paek et al., 2019). During summer, the position of WPSH shifts eastward and weakens in intensity during El Niño events, whereas during La Niña events, the WPSH position shifts westward and its intensity strengthens (Liu et al., 2023). This phenomenon suggests that the air-sea interaction under the influence of the subtropical high may play a crucial role in SST rise (Zeng et al., 2010). Additionally, it may be a key factor contributing to higher SST during summers characterized by La Niña events compared to those with El Niño events (Timmermann et al., 2018). In the northern part, possibly due to significant disturbances from local factors (such as Q_LH), although MEI and T_max exhibit a negative correlation, the correlation is not significant ( Table 1 ).

Based on the interannual variation trend of SST, the causes of extreme high-temperature events in certain years need to be analyzed specifically (Gao et al., 2020; Yan et al., 2020; Tan et al., 2023). In both the southern and northern parts of the Eastern China Seas, T_max in 2016 was the highest in 38 years ( Figure 5 ). Several studies have reported extreme high-temperature events in the Eastern China Seas during the summer of 2016 (Tan and Cai, 2018; Moon et al., 2019; Yan et al., 2020). In August 2016, the region experienced an exceptional warming event with maximum anomalies exceeding 2°C (Tan and Cai, 2018), and the intense warming led to SST breaking historical records observed by satellites (Yan et al., 2020). The high SST in 2016 is attributed to various factors such as increased net heat flux and enhanced warm advection (Tan and Cai, 2018), increased Changjiang diluted water (Moon et al., 2019), and a high-pressure system over the Eastern China Seas splitting from the western Pacific subtropical high (Yan et al., 2020). Gao et al. (2020) attributed the high temperature events during the consecutive summers from 2016 to 2018 to shortwave radiation and ocean advection anomalies, as well as the reduction of surface wind speeds. The high-temperature event in 2022 was caused by a combination of anomalous atmospheric and oceanic conditions, including enhanced insolation, weakened surface wind speed, suppressed latent heat loss from the ocean, a shallower mixed layer, and an upper ocean current anomaly (Tan et al., 2023).

Understanding the predictability of T_max is crucial for assessing the impacts of extreme heat events on marine ecosystems and coastal communities. In the southern part of the Eastern China Seas, T_max is closely related to AirT_max and its timing, suggesting that AirT_max can serve as a reliable predictor for T_max. This relationship is more complex in the northern part, where Q_LH significantly influences T_max timing. High Q_LH during August tends to suppress further increases in SST, advancing T_max timing. The negative correlation between Q_LH and T_max timing (-0.56, p<0.01) indicates that Q_LH can be used to correct T_max timing forecasts. Therefore, a hybrid approach integrating AirT_max and Q_LH can enhance the accuracy and timeliness of T_max forecasts in the northern part. Furthermore, other parameters influencing SST, such as ocean currents, wind, solar radiation flux, MEI, and so on, can be integrated to further improve the predictive capability. Machine learning algorithms, can handle the non-linear relationships and interactions among multiple parameters (Dong et al., 2021; Manucharyan et al., 2021), thereby providing robust predictions of T_max. The impacts of these trends on marine ecosystems and fisheries can be profound. For example, higher SSTs can lead to shifts in species distribution, affecting local fisheries (Hobday et al., 2016; Oliver et al., 2018). Specific case studies, such as the 2016 extreme high-temperature event, highlight the importance of understanding these dynamics (Tan and Cai, 2018; Moon et al., 2019; Yan et al., 2020). Integrating various parameters and employing machine learning techniques can further enhance the precision and lead time of T_max forecasts, facilitating better preparedness for extreme heat events and their impacts on marine environments. This is a topic that is worth studying further.

While our study provides significant insights into the patterns and trends of T_max in the Eastern China Seas, several limitations must be acknowledged to offer a balanced perspective on our findings. Our initial analysis indicates that the faster warming rate of T_max in the northern region is largely due to the corresponding trends in AirT_max, which also exhibit a faster warming rate in the north. Additionally, the differing trends in T_max timing between the northern and southern regions are closely linked to the distinct trends in Q_LH during July and August. However, the underlying mechanisms behind these trends remain complex and are not fully explored in this study. Future research should delve deeper into understanding the reasons behind the varying trends in AirT_max and Q_LH between the northern and southern regions. By addressing these complexities, future studies can build upon our findings to provide a more comprehensive understanding of the regional climate dynamics affecting T_max. Moreover, while our study focused on the observed trends and initial causal factors, further investigation into the broader climatic phenomena and their interactions with local processes will be crucial in refining our understanding of T_max variability in the Eastern China Seas.

Our study reveals a significant positive correlation between T_max and the probability of MHWs occurrence in the Eastern China Seas ( Figure 8 ). The average probability of MHWs occurrence in the study area is approximately 6%. However, during T_max, this probability surges to 16%, indicating a heightened risk of MHWs coinciding with T_max. The year 2016 exemplifies this relationship, as it recorded the highest T_max and the most frequent MHW occurrences ( Figure 8 ). In 2016, MHWs affected 69% of the study area during T_max, meaning that 69% of the region experienced MHWs when SST reached its annual maximum. This dramatic increase underscores the strong linkage between T_max and MHWs.

The implications of this relationship between T_max and MHWs are profound. MHWs are known to have devastating impacts on marine ecosystems, including widespread coral bleaching, altered species distributions, and disrupted marine food webs (Hobday et al., 2016; Oliver et al., 2018). The concurrent occurrence of T_max and MHWs amplifies these adverse effects, posing severe threats to marine biodiversity and the health of coastal environments. Recognizing the critical role of T_max in triggering MHWs, it is essential to monitor and predict T_max events accurately. Enhanced forecasting of T_max can provide early warnings for impending MHWs, enabling timely mitigation efforts to protect vulnerable marine habitats. Additionally, understanding the mechanisms linking T_max with MHWs can inform strategies to enhance the resilience of marine ecosystems to extreme thermal events.