Three sources of satellite data were employed in this study, including the Moderate Resolution Imaging Spectroradiometer (MODIS) -Aqua Daytime SST, Multiscale Ultrahigh Resolution (MUR) SST and Himawari-8 SST. The available periods of these data are 2002-07 to ongoing, 2002-06 to ongoing and 2015-06 to ongoing, respectively.

Monthly averaged MODIS -Aqua SST Level 3 Mapped data with a spatial resolution of 4 km × 4 km during 2003 to 2021 were obtained from the U.S. National Aeronautics and Space Administration (NASA) website ¹ . Together with MODIS-Aqua SST, two more datasets were used for comparison and verification. The MURSST Level 4 SST, produced by NASA JPL Physical Oceanography DAAC ² with a 0.01-degree resolution, is based upon nighttime GHRSST L2P skin and subskin SST observations from several instruments, including the AMSR-E, AMSR-2, MODIS, WINDSAT radiometer, AVHRR, and in situ SST observations from the NOAA iQuam project. In this study, daily MURSST data from October 2020 to June 2021 are used. The monthly Level 3 Himawari-8 SST data with a 2 km spatial resolution from October 2020 to June 2021 were provided by the Japan Meteorological Agency (JMA) ³ , this product is mainly based on the observations by the new generation geostationary satellite Himawari-8. Notably, the Himawari-8 SST has been validated against the in-situ SST data (iQuam) and shows reliable quality.

The seafloor observations are achieved through the submarine cable observation systems (Zhai et al., 2020). These systems were applied to online monitor temporal variations in marine environment of the ecologically important regions. Four submarine cable observation systems are marked off with red squares around the Changjiang and Qiantangjiang estuaries in Figure 1B . The submarine cable online observation systems were placed at 20 cm above the sea floor and equipped with a HydroCAT conductivity-temperature-depth (CTD) recorder to continuously observe temperature, salinity, pressure and other variables in the bottom layer. In this study, hourly averaged temperature in bottom layer and the water elevation were used to examine the performance of the numerical model.

ERA5 ⁴ is the fifth generation ECMWF reanalysis of the global climate and weather, which provides hourly estimates for a large number of atmospheric, ocean-wave and land-surface quantities. In this study, 0.25-degree hourly wind vectors and surface radiation flux data, including Surface net solar radiation (SNSR), Surface net thermal radiation (SNTR), Surface sensible heat flux (SSHF), and Surface latent heat flux (SLHF), of 2020 to 2021 are used as the atmospheric forcing in our model. Meanwhile, ERA5 data from 2003 to 2021 were also utilized in the correlation analysis in section 4. The ECMWF convention for vertical fluxes is positive downwards, and net heat flux from ERA5 is calculated following:

SST gradient magnitude was calculated by the central difference method as:

Where ∇ _x SST and ∇ _y SST are the zonal and meridional components of the SST gradient, ΔX and ΔY are the distances in kilometers between the neighboring grid points in the longitudinal and latitudinal directions, respectively. |∇SST| is the gradient magnitude. This algorithm has been widely used in the calculation of coast SST gradient (Vazquez-Cuervo et al., 2013; Saldias and Lara, 2020; Saldias et al., 2021).

The studied region characterizes irregular coastline and rugged topography, therefore the traditional circulation model working under structured grid is no longer capable. For this reason, the Semi-implicit Cross-scale Hydroscience Integrated System Model (SCHISM) was employed in this study (Zhang et al., 2016a; Zhang et al., 2016b). SCHISM is a 3D unstructured-grid model developed from the original Eulerian-Lagrangian Finite Element model, using the semi-implicit time stepping and hybrid finite-element/finite-volume method with the Eulerian-Lagrangian algorithm to solve the Navier-Stokes equations in hydrostatic form. In addition, it provided a new hybrid SZ coordinate called Localized Sigma Coordinates with Shaved Cells (LSC²) in the vertical dimension, which makes no bathymetry smoothing possible. More details about the SCHISM model can be found on its official website ⁵ .

The model region covers the whole Bohai sea, Yellow sea, and East China Sea with three open boundaries ( Figure 1A ). In order to better simulate the hydrodynamics in the nearshore areas and save computational resources, the locally-refined grid was designed. The horizontal resolution near the open boundaries is depth-dependent varying between 7 to 50 km, while for the refined regions (blue dotted box in Figure 1A ), the horizontal resolution ranges from 2 km in the shallow coastal water to 8 km in the deeper area. The whole mesh grid contains 24742 nodes and 47319 triangle elements. In the vertical direction, 59 layers LSC² grid was applied and the bathymetry was interpolated from the Shuttle Radar Topography Mission (SRTM) 15Plus V2.1 (Tozer et al., 2019).

The 3-D water temperature, salinity, velocity, and sea surface height for initial and boundary conditions were interpolated from the 41-layer Hybrid Coordinate Ocean Model and Navy Coupled Ocean Data Assimilation (HYCOM + NCODA) Global 1/12° Analysis data. Additionally, we prescribed the tidal forcing of eight tidal components, including M₂, S₂, K₂, N₂, K₁, O₁, Q₁, and P₁ at the open boundaries. All the tidal information comes from the latest version of the Finite Element Solution (FES2014) tide model by AVISO+ ⁶ . We considered two major river sources, the Yangtze River and the Qiantang River, the monthly averaged river discharge data is from the River and Sediment Bulletin of China ⁷ , besides, the temperature of the rivers is based on MODIS-Aqua 4km daily SST product, and the salinity is set to 3 PSU uniformly. Hourly atmospheric forcing was generated using the ERA5 product originated from Climate Data Store⁴. The simulation period was initialized on October 1, 2020 and ended on December 31, 2021. During the integration, the time step was 120s and hourly outputs were saved.

Figure 2 displays the monthly mean SST gradient magnitude and the 10m-high wind vectors from October, 2020 to March, 2021 within the study area. In November, the ZFCF started to form under the strong prevailing northeast monsoon, then its intensity kept strengthening along the 50-m isobath and peaked in January, with a maximum intensity of 0.2°C/km. In the following February and March, there was a significant decrease in both the intensity and spatial range of the ZFCF, especially in the area south of the 28°N, the SST gradient has dropped by more than 0.12°C/km and the front almost disappeared by March. However, compared with the evolution of winds in January, there is no significant change in the local wind both in speed and direction.

To show this phenomenon in more detail, we calculated the time series of spatially-averaged daily SST gradient as well as the local wind vectors of the focus region (blue dashed box in Figure 1A ) and presented them in Figure 3 . As revealed by Figure 3 , when entering January, the SST gradient shows a high correlation with the southward meridional wind, e.g., it increases significantly during these two strong northerly wind processes around early January and mid-February, respectively, which agrees with the points of previous works (Cao et al., 2021). Then an abnormal abrupt decrease of the SST gradient occurred: the magnitude of SST gradient in the focus region quickly dropped by more than 0.05 °C/km during 15 days after late January, which is highlighted in red in Figure 3 . This phenomenon can be divided into two periods: in the first period, the northerly wind significantly weakens and there is a synchronous dramatic decrease in SST gradient; in the second period, after the recovery of the northerly monsoon, the SST gradient does not recover as expected but oppositely a repaid decline process occurs. So, how does this phenomenon develop?

The satellite observed SST has captured the rapid changes of the ZFCF successfully, however, it is hard to further investigate the regulatory mechanism of this phenomenon, thus the numerical model is employed to make further explorations. To ensure the reliability of the discussions, we first validated the model performance in describing local hydrodynamics.

The vertical structure of temperature and current across the ZFCF in two representative sections (sec1 and sec3, Figure 1B .) from October 2020 to March 2021 are shown in Figures 8 and 9 , respectively. As Figures 8A-F reveal, in October, despite the presence of strong ZFCC, no significant thermal fronts occurred because the water column is well stratified. As it comes to November, the whole water column is sufficiently mixed under the strong wind stress, meanwhile the cold ZFCC dominates the coastal area, therefore the horizontal temperature gradient begins to strengthen and an apparent thermal front (ZFCF) forms at regions 20~50 km off the coast. After reaching its peak in January, the vertical isotherm of the ZFCF gradually shifts direction and transforms into a horizontal distribution, suggesting the decay of ZFCF. From Figures 8G-L , the boundary of the northward and southward current can be clearly distinguished by the black dotted lines; interestingly, the boundary of two reversed current is located about 60 km off the coast, which is remarkably different with the location of ZFCF. This feature seems to be contradictory to the formation mechanism of the front raised by former studies, which will be explained in section 4. The first branch of KI carries relatively warm water and joins into the ZFCC near the coast. Thus, a front closer to the coast is formed. Meanwhile, the weakening of northward flow in January also suggests that the strong KI and northly monsoon could inhibit the southward flow within the Taiwan Strait. Obviously, Figure 8 provides a piece of clear evidence that the ZFCF is affected by the KI.

Sec 3 is located at the junction of the Zhejiang and Fujian coasts. Along sec 3 ( Figures 9A-F ), the temperature structure of the upper water has the same pattern as sec 1 in regions shallower than 50 m, while the deep water has strong stratification during October to December. The water column could be well mixed only in January, and then it stratifies again accompanied by the weakening of the front. The variation in the stratification of the water column is usually dominated by the heat fluxes, and thus exerts potential impact on the ZFCF, which will be discussed in the next section. In addition, the flow structure reveals that the southward (ZFCC) and northward (TWC) flows have a clear decreasing tendency through February and March, which is a strong indicator that the change of background flow may be the direct cause of the abrupt change of ZFCF.

The intensity of the ZFCF is closely related to the northeast monsoon ( Figure 2 ): the monsoon weakened to some extent during the sudden drop of the ZFCF. Further, the stratification of the water column caused by variation of the net heat flux may also contribute to the sudden drop in thermal front ( Figures 8 and 9 ). Therefore, based on correlation analysis, the relationships between the heat flux, wind stress, and the ZFCF were analyzed. We calculated the area-averaged local winds, net heat flux, and the SST gradient from January to March through 2003-2021 and the results are presented in Figure 10 . As displayed, the intense heat loss (negative net heat flux) corresponds well with higher SST gradients, the regression analysis result (R=-0.72, P<0.001) also confirms their significant correlation. During a heat loss process, the shallow nearshore area tends to be cooled down more completely compared with the deep offshore region due to stronger vertical mixing and relatively lower heat content over the whole water column. Therefore, the temperature difference between nearshore and offshore waters gets further expanded, triggering a stronger ZFCF.

Winds can directly modulate the circulation, vertical mixing, and the air-sea heat exchange process (Qiao et al., 2006; Wang et al., 2015; Zhou et al., 2015; Park et al., 2020), further adjusting the ZFCF. Thus, as speculated, Figure 10 shows significant correlationship between SST gradient and wind stress (R=0.75, P<0.01). However, effect of winds on ZFCF is a comprehensive process and next we will explore the contribution of winds from different aspects.

The encounter of two different water masses or currents is one of the major frontogenesis mechanisms, as section 3.3 reveals, the weakening of ZFCC and TWC is highly correlated with the abrupt decay of ZFCF. We used the monthly-averaged HYCOM current product and MODIS-Aqua SST data to find out the relationship between the currents and ZFCF from January to March over an interannual scale. Expect for sec 1 and 3 presented above, two more sections (sec 2 and 4, Figure 1B ) are selected to investigate the role of KI. The heat flux Q_T through four representative sections are calculated to diagnose the interactions between the background flow and ZFCF following:

where T is the temperature of the water column, V is the cross-section velocity, η is sea water level, h is the water depth, and l is the width of the section. The southward current through sec 1 and 3 represent the ZFCC, while the northward current across sec 1 represents the flow through Taiwan Strait, and the westward currents crossing sec 2 and 4 represent the KI.

Figure 11 shows the time series of area-averaged SST gradient and Q_T through sec 1-4, meanwhile, results of linear regression analysis (at 95% confidence level) are shown in Table 3 . For the ZFCC (southward current through sec 1 and 3), the variation of Q_T and SST gradient synchronize well with each other and significantly correlated (R=0.46 and 0.56, respectively). The ZFCC flows southward along the coast, bringing cold and fresh water to the focus area, and therefore directly maintains and regulates the ZFCF.

For the warm flow through Taiwan strait (northward current across sec 1), it leads to a significantly negative correlation between the Q_T and SST gradient (R=-0.66, P<0.001). This seems to be inconsistent with common sense, as the encounter of two different flows with larger temperature gradients tends to establish a strong thermal front. And here, the Taiwan Strait current and ZFCC have opposite directions, strong northerly winds will strengthen southward ZFCC and weaken northward Taiwan Strait current, thus the Q_T of Taiwan Strait current through sec 1 is reducing. The ZFCC and TWC could be regarded as the cool and warm water sources for forming ZFCF, when TWC weakens, extra warm water sources are needed to maintain the strength of ZFCF. For the ZFCF, KI is the possible one which might supply warm water except for the Taiwan Strait current.

Wang and Oey (2014) indicated that the strengthened northeast wind would significantly enhance the KI. Recently, Kang and Na (2022) further reported significant correlation relationships between KI and wind stress and wind stress curl. Does this mean stronger northeast winds will generate stronger KI, and thus leads to a stronger ZFCF? The answer is uncertain because stronger northeast winds enhance the KI but will inhibit the northward transport of warm KC water. This is why Figure 11 and Table 3 both suggest no significant relationship between the KI (westward current in sec 2 and 4) and the SST gradient.

Discussions above point out that the ZFCF is closely modulated by the heat flux and background circulation in long term. To find out the key factors leading to this abrupt change of ZFCF in early spring, 2021, heat budget analysis was conducted based on numerical results. Since the onshore water column was well mixed during the abrupt change, we selected a series of representative points (blue and red dots in Figure 12E ) located on either side of ZFCF, and the average heat budget through the whole water column at each point was calculated following Vijith et al. (2020) and the equation was written as:

Where T, ρ₀, c_p and h represent the temperature, mean density, water column depth and specific heat capacity of seawater, respectively. u,v and w are the zonal, meridional and vertical velocity. The suffix a are the vertically-averaged quantity and the subscript -h denotes the quantity at the base of the water column. q₀ and q_pen are the surface net heat flux and penetrative loss of shortwave radiation, respectively. K_H and K_z represent the horizontal and vertical eddy diffusivities, respectively.

The sum of the horizontal advection, horizontal mixing, vertical mixing, entrainment and net heat flux well consist the overall pattern of the tendency term with a correlation of 0.87 and 0.83, RMSE of 0.13 and 0.09 °C/day (within 95% confidence interval) in cold water and warm water, respectively. The result is a satisfyingly good closure of the whole water column heat budget for such short time scales. Our analysis points out that net heat flux is important in determining the net tendency of the temperature on both sides of the front ( Figures 12A and B ), while in the cold water west to the front, vertical mixing and horizontal advection also greatly contribute to the tendency before February. For the abrupt change around later January, the temperature tendency in cold water has already been positive, while that in warm water east to the front was still negative. That indicates the temperature in cold water begins to warm up thus the temperature difference between two water masses is reduced. When the positive tendency reached its peak (~0.18°C/day) in cold water, the intensity of the ZFCF starts abruptly decreasing, and the net heat flux was the predominant contributor.

Interestingly, during the change period, the contributor of the positive tendency changed from net heat flux to horizontal advection in cold water, which provides a good answer to the question at the end of section 3.1. Although the recovery of the monsoon causes the net heat flux tends to cool the onshore water again, the heat carried by the horizontal advection compensated for the loss and kept the tendency positive, then resulting the decrease of ZFCF in the second period. Besides, the decomposition of the horizontal advection reveals that the flow pattern might change during the abrupt change of the ZFCF, because the zonal advection term changed from positive to negative and the meridional advection changed from negative to positive on both sides of the front. Combined with the analysis of the above sections, the variations of the horizontal advection of ZFCC and KI might induce the positive temperature tendency.

To further illustrate the contributions of different dynamic factors in triggering the abrupt decay of ZFCF in February, we conducted a series of sensitive numerical experiments based on above analyses. Each experiment runs for seven months from October 1, 2020, the control conditions of these sensitive experiments are listed in Table 4 .

The differences in SST gradient between the sensitive experiments and control run are shown in Figure 13 . When the wind stress forcing is closed (No-wind case), the SST gradient have changed dramatically in both distribution and intensity, the ZFCF almost disappears; in this case, thermal fronts mainly occurred around the Yangtze River estuary or regions with complicated bathymetry as a result of river discharge and modulation of current by topography. Moreover, the fronts triggered by the KI to the northeast of Taiwan Island are further away from the shelf, indicating a decreased KI under weak northerly monsoon, which agrees well with previous study (Wang and Oey, 2014). We found that the ZFCC almost disappeared with a small velocity lower than 0.1 m/s; thus, local wind stress is a key factor during the formation of ZFCF through modulating the intensity of ZFCC. When the wind stress in February and March were replaced by that in January (Exp. Wind-Jan, Figure 13E ), the ZFCF was significantly enhanced, again proving that the weakening of the northeast wind stress plays a vital role during the abrupt decay process in February and March.

When no river discharge was considered (Exp. No-river, Figure 13B ), the location of ZFCF is closer to the coastline in January compared with the control run, besides, the intensity of ZFCF decreases slightly in February and March. As an essential component of the ZFCC, the riverine flow merges with the ZFCC, increasing the flow flux and strengthening the zonal width of the ZFCC (Wu et al., 2021), which then leads to a farther ZFCF off the coast. Besides, buoyancy forcing produced by the Yangtze River discharge can strengthen the southward flow of the Zhejiang-Fujian coastal water (Wu et al., 2021; Zhang et al., 2022), so the ZFCF decreases in February and March under the condition of no river discharge. Further, runoff of the Yangtze River reduced from January (40.6×10⁹ m³) to February (33.53×10⁹ m³), which might also contribute to the decrease of the ZFCF in February.

From experiments No-tide and Rad-Jan ( Figures 13C and D ), limited effect of solar radiation and tidal forcing on the intensity of ZFCF could be speculated, and the front is fixed on the topographic slope under the tidal mixing (Cheng et al., 2017). When tidal forcing is neglected, the location of the front moves onshore. Moreover, the result also points out that the net heat flux, which is important in determining the net temperature tendency, is mainly modulated by latent and sensible heat flux controlled by winds.

To investigate the influence of the KI, Exp. No-KI was designed by blocking the western passage of the KI ( Figure 13F ) to simulate the situation of the weak KI. Notably, KI and the shelf circulation are vital components of the ECS, to balance the mass flux in and out of the continental shelf, when the KI was cut off, the currents through the Taiwan Strait will be influenced (Yang et al., 2020; Wu, 2021). In this case, the front formed by the KI in northern Taiwan Island has disappeared completely, meanwhile, the ZFCF has significantly changed. The reduction of frontal intensity in January is much weaker than that in February and March, this phenomenon indicates that the KI can hardly reach the focus region under the strong northeast monsoon in January, so the reduction of ZFCF in the southern area is much remarkable than that in the northern area. Tracer experiments by Yang et al. (2018) suggested that the Kuroshio subsurface water cannot reach the algal bloom area along the Zhejiang coast in winter, when the prevailing monsoon relaxes, the onshore intrusion of the Kuroshio will be less prohibited by the barotropic pressure force created by the onshore Ekman transport (Yang et al., 2012; Xu et al., 2018). The intensity of ZFCF decreases distinctly in February under the absence of the KI, the maximum reduction reaches 0.08 °C/km inside the focus region, which is comparable with the change in the No-wind experiment, and even larger in March. Exp. No-KI indicates that the ZFCF is formed by the southward advection of the ZFCC, which was mainly forced by the northeast wind forcing, and the TWC water. As Figure 13F revealed, only a minor change of ZFCF occurred in January when the KI is cut off. However, this circumstance changed dramatically in Feb and Mar when the prevailing monsoon relaxed. This is attributed to the impact of KI, which can reach the focus region and then further affect the ZFCF more efficiently under the weak northeast monsoon. At this point, assisted by the heat budget analysis, the mechanism of the decrease in the second period can be revealed. During the weak of ZFCF in the first period which was caused by relaxation of the monsoon, the invaded water tended to transport to the focus region at the same time. And after this period of monsoon relaxation, the invaded water reached the focus area and heated the weak ZFCC water, so the component of horizontal advection term contributed positive temperature tendency changed from zonal shelf warm water advection to meridional KI advection ( Figures 12C and D ). Although, the monsoon recovered and tended to cool the onshore water again, the heating effect by KI was stronger than wind cooling effect, thus ZFCF was continually weakening in the second period until the winds forced the cool ZFCC southward and impeded the northward transport of the KI water again. However, Figure 3 shows that the intensity of ZFCF gets stronger after the abrupt decrease, which is caused by the interaction of cool water forced by recovered winds and the residual warm KI water left on the shelf. Also, in the following period with relatively weak northeast monsoon, the continuously invaded water maintained the less intense ZFCF until the onshore water was heating up in the middle March.