In this study, we used the monthly average SST measured by MODIS-Aqua (https://oceandata.sci.gsfc.nasa.gov/directaccess/MODIS-Aqua/Mapped/Monthly/4km/sst/) from 2003–2021 with a spatial resolution of 4 km. Low-quality data were eliminated according to the quality control policy of the MODIS products (quality levels>0). Meanwhile, due to the apparent distortion of the SST in the shallow water region of the nearshore bay, we also eliminated the SST where the water depth was less than 2 m using the ETOPO1 topography (https://sos.noaa.gov/catalog/datasets/etopo1-topography-and-bathymetry/). The SST data from June to August were averaged to obtain the summer SST data.

Since the SST variation is influenced by atmospheric conditions, the atmospheric reanalysis data is analyzed in this study. The monthly sea surface forcing (including 10 m wind speed, air-sea heat flux and 2 m air temperature) from 2003–2019 were derived from the ERA-interim reanalysis dataset (https://apps.ecmwf.int/datasets/data/interim-full-moda/levtype=sfc/) with a horizontal resolution of 0.125°. Monthly precipitation data were obtained from the Integrated Multi-satellite Retrievals for GPM (IMERG) (https://disc.gsfc.nasa.gov/datasets/GPM_3IMERGM_06/summary), which had a horizontal resolution of 0.1° (Huffman et al., 2019). When exploring the relationships between SST variation and sea surface forcing, the corresponding local forcing are used.

To analyze the influence of the sea surface wind on upwelling, the wind stress and wind stress curl were calculated. The wind stress was calculated using the quadratic law:

where, τ → =τ _x i+τ _y j is the wind stress vector, τ _x and τ _y are the zonal and meridional wind stresses, respectively. ρ_a is the air density, V → = u i + v j is the 10 m wind vector, u and v are the zonal and meridional components of the wind velocity, respectively. C_d is the drag coefficient whose calculation formula is as follows (Yelland and Taylor, 1996; Yelland et al., 1998):

Since Ekman transport is a key dynamic process inducing the Qiongdong upwelling, we also calculated the along-coast wind stress component in the eastern part of Hainan Island (the angle between the coastline and the north direction is approximately 30°).

The wind stress curl is calculated using the wind stress:

Empirical orthogonal function (EOF), also known as eigenvector analysis, is a method for analyzing the structural features in matrix data and extracting the main data features. EOF has been widely used to analyze temporal and spatial variations in atmospheric and oceanic space fields (Dong et al., 2018; Syamsuddin, 2020). In this study, we used the EOF to obtain the main modes of temporal and spatial variations in SST around Hainan Island (18–21.5°N, 108–112°E) and analyze its interannual variation characteristics. Univariate linear regression was used to fit the long SST trend in the upwelling region. Interannual SST anomalies were obtained by subtracting the SST regression values from the corresponding SST observations. The Pearson correlation analysis method was used to calculate the correlation coefficient and its significance level when analyzing the correlation between different factors.

Upwelling along the coast of Hainan Island occurred almost every summer. The average summer SST from 2003–2021 is shown in Figure 1B . The Qiongdong upwelling is located on the east side and northeast side of Hainan Island (18–21°N, 110–112°E). The SST over the Qiongdong upwelling region was 1–2°C lower than the surrounding waters. The southernmost part of the low-temperature region of Qiongdong upwelling reaches the coastal waters of Sanya City and the northernmost part reaches a latitude of 21°N latitude. The Qiongxi upwelling on the west side of Hainan Island has a relatively small area (18–19.5°N, 108–109°E) with a SST of 0.5–1°C lower than the surrounding waters. The Qiongdong upwelling covers a larger area than the Qiongxi upwelling, and the SST over the Qiongdong upwelling region is lower than that over the Qiongxi upwelling region. This indicates that the Qiongdong upwelling was stronger than the Qiongxi upwelling.

The SST in Qiongdong and Qiongxi upwelling regions have evident interannual variations. To explore SST changes in the two upwelling, the region surrounded by the 29.0°C isotherm on the east side of Hainan Island (A_east ) is used to represent the waters affected by the Qiongdong upwelling, and the region surrounded by the 29.5°C isotherm on the west side of Hainan Island (A_west ) represents the waters influenced by Qiongxi upwelling ( Figure 1B ). The two areas (A_east and A_west ) and the mean temperature within the marked areas ( S S T ¯ e a s t and S S T ¯ w e s t ) were calculated to characterize surface low-temperature water.

Figure 2 shows the interannual variations and long-term trends of SSTs in Qiongdong and Qiongxi upwelling regions during 2003–2021. This shows that both the mean temperature and area affected by upwelling showed evident interannual variations. For the Qiongxi upwelling, the S S T ¯ w e s t showed a significant increasing trend of 0.012°C/year (p< 0.01) and the A_{wes t} showed a significant decreasing trend of -428 km²/year (p< 0.05) ( Figure 2B ). The interannual changes in S S T ¯ w e s t and A_west were significantly correlated (r = -0.80, p< 0.01). As for the Qiongdong upwelling, the S S T ¯ e a s t showed an increasing trend of 0.0047 °C/year (p = 0.30) and the A_east showed a decreasing trend of -401 km²/year (p = 0.28) ( Figure 2A ), although the trend was not significant. In addition, the interannual variation in S S T ¯ e a s t and A_east showed a correlation coefficient of -0.74 (p< 0.01). Therefore, the SSTs in both the Qiongdong and Qiongxi upwelling regions exhibit warming trends, and the tendency is more pronounced in the Qiongxi upwelling region than in the Qiongdong upwelling region, as SST trend in Qiongxi upwelling region is higher than that in Qiongdong upwelling region. Meanwhile, the areas of the surface low-temperature regions exhibit decreasing trends.

Correlation analysis was conducted to compare the SST in Qiongdong and Qiongxi upwelling regions. It was found that the A_east and A_west were large and fluctuated significantly before 2013, whereas they were small and relatively stable after 2013 ( Figure 2 ). The interannual changes in A_east and A_west were significantly correlated (r = 0.80, p< 0.01). However, there was no correlation between the interannual changes in S S T ¯ e a s t and S S T ¯ w e s t (r = 0.22, p = 0.36).

For non-upwelling sea regions, we calculated the average SST in surrounding sea region E ( S S T ¯ r e g E ) and region W ( S S T ¯ r e g W ) which were remarked in Figure 1B . Results show that S S T ¯ r e g E and S S T ¯ r e g W have increasing trends of 0.036°C/year (p< 0.05) and 0.055°C/year (p< 0.01), respectively. It suggests that the SST warming trend in the upwelling regions is smaller than that in the surrounding sea region ( Figures 2A, B ).

As shown in Figure 1B , the Qiongdong upwelling has a large range and can be further divided into two parts (east and the northeast) by a boundary at a latitude of 19.64°N around Wenchang City. In this subsection, we used the lowest SST in the upwelling region to identify the upwelling core. The location and water temperature of the upwelling cores on the east (ES), northeast (NES), and west (WS) sides of Hainan Island were calculated ( Figure 1B and Figure 3 ). The locations of the three cores were relatively stable. The ES upwelling core was at 19.01°N, 110.58°E; the NES upwelling core was at 20.12°N, 110.95°E; and the WS upwelling core was at 18.62°N, 108.54°E.

As for the long-term trend, the latitudes of the ES and WS upwelling cores fluctuated slightly. However, the core of the NES upwelling tended to move south ( Figure 3A ). Figure 3B shows the interannual variation in water temperature of the WS, ES, and NES upwelling cores. The mean temperature of the WS upwelling core was 28.68 °C with a warming trend of 0.047 °C/year (p< 0.01). However, for the Qiongdong upwelling, the mean temperature of the NES upwelling core was 27.80°C with a warming trend of 0.035 °C/year (p = 0.06), and the mean water temperature of the ES upwelling core was 27.60°C with no significant long-term trend.

We also selected three main upwelling regions (red boxes in Figure 1B ) to study the statistics of interannual SST variation. As shown in Figure 4A , the SSTs in the NES and WS upwellings show significant warming trends of 0.036°C/year (p< 0.05) and 0.045°C/year (p< 0.01), respectively. The mean SST in the ES upwelling showed a warming trend of 0.010°C/year, although this was not significant (p = 0.61). Notably, the SST warming trends of the main upwelling regions were similar to those of the upwelling cores ( Figure 3B ).

In addition to the long-term trends, SST has an interannual variation. we calculated the detrended SST anomaly (SST minus its corresponding trend, dSST anomaly) to represent the interannual variation over the ES, NES, and WS upwellings. As shown in Figure 4B , the dSST anomaly in three upwelling regions show evident interannual variation, which has an amplitude of 1.57°C, 1.33°C, and 0.94°C in the ES, NES, and WS upwelling regions, respectively. After removing the warming trend ( Figure 4A ), the interannual variation in the dSST anomaly around Hainan Island remained consistent. Specifically, the correlation coefficient between the ES and NES upwelling regions is 0.60 (p< 0.01), between the ES and WS upwelling regions is 0.52 (p< 0.05), and between the NES and WS upwelling regions is 0.63 (p< 0.01). Although the dSST anomalies in the three upwelling regions were correlated with each other, the dSST anomalies in the ES and NES upwelling regions still had some differences. The ES upwelling had lower dSST anomalies in 2012 and 2018; however, the same cold anomalies did not appear in the NES upwelling ( Figure 4B ). The relationship between the ES and NES upwelling regions is further explored in the following section.

To further analyze the temporal and spatial characteristics of the Qiongdong and Qiongxi upwellings, we carried out an EOF analysis to obtain SST’s spatial mode and the associated time series (PC) of SST. Figure 5 shows the spatial distribution and PC values of the first three modes. The variance contribution rates of the first three EOF modes are 52.57% (EOF1), 10.47% (EOF2), and 5.87% (EOF3). Correlation analysis was performed for PC1, PC2, and PC3 with SST and dSST anomalies in the ES, NES, and WS upwelling regions. The results are presented in Table 1 .

The first EOF mode had a positive SST pattern and exhibited a consistent warming trend ( Figures 5A, B ). To illustrate the connection with global warming, the global average temperature and SST anomalies (GLBTSSST) index (https://psl.noaa.gov/gcos_wgsp/Timeseries/) provided by NOAA was used. The PC1 and GLBTSSST indices had a correlation coefficient of 0.77 (p< 0.01), indicating that PC1 represents global warming. Meanwhile, PC1 had a significant positive correlation with the interannual SST variation in the ES, NES, and ES upwelling regions ( Table 1 ), which reflects the warming trend of the three upwelling regions ( Figures 2 , 3B and Figure 4A ). The increase in PC1 is consistent with the area reduction of the upwelling regions after 2013 ( Figures 2A, B ), suggesting that global warming promoted the shrinkage of low-temperature region affected by upwelling around Hainan Island. According to the correlation coefficient between PC1 and the SSTs in the three upwelling regions ( Table 1 ), SST in WS upwelling region has a stronger response (r = 0.85) to global warming than ES upwelling (r = 0.63) and NES upwelling (r = 0.67). These results are consistent with the conclusions of the previous section ( Figures 3B and 4A ).

The second and third modes show the spatial heterogeneity of the SST. It was positive over the NES upwelling region for the second mode ( Figure 5C ), while it was positive over the ES upwelling region and part of the NES upwelling region for the third mode ( Figure 5E ). Correlation analysis showed that PC2 ( Figure 5D ) was significantly correlated with the interannual variations in SST (r = 0.51, p< 0.05) and dSST anomaly (r = 0.71, p< 0.01) over the NES upwelling region. In contrast, PC3 ( Figure 5F ) is significantly correlated with the interannual variations of SST and dSST anomaly over the ES upwelling regions with correlation coefficients of 0.50 (p< 0.05) and 0.55 (p< 0.05), respectively. Lin et al. (2016b) suggested three types of upwelling patterns: intensified ES upwelling, intensified NES upwelling, and both the ES and NES upwelling. We found that the intensified NES upwelling pattern was similar to the spatial distribution of EOF2 ( Figure 5C ); meanwhile, the intensified ES upwelling pattern was similar to that of EOF3 ( Figure 5E ). For example, the Qiongdong upwelling in 2015, 2020, and 2021 showed a similar spatial distribution to EOF2 with a weaker ES upwelling but a stronger NES upwelling. However, in 2012 and 2018, the Qiongdong upwelling had stronger ES upwelling but weaker NES upwelling, showing a similar spatial distribution to that of EOF3 (not shown). Therefore, we can conclude that the second EOF mode mainly reflected the interannual SST variation in the NES upwelling region, while the third EOF mode mainly reflected that in the ES upwelling region. The independence between EOF2 and EOF3 reveals that different factors influence the interannual SST variations in the ES and NES upwelling regions.

Low SST is always used to identify the occurrence of upwelling. However, it cannot directly reflect the strength of upwelling as it is also influenced by atmospheric condition. From satellite-derived SST, both the SST in the Qiongdong upwelling and Qiongxi upwelling regions show a warming trend during 2003-2021. Under the background of global warming, a warming SST does not mean a weakened upwelling. SSTs in Chinese marginal seas have been widely proved warming (Li et al., 2017; Li et al., 2019; Liu et al., 2020), so is the SST in the northwest area of the South China Sea ( Figure 5A ). To explore the change trend of upwelling, we compared the SST between upwelling region and surrounding non-upwelling region. As SST warming is evident in both upwelling regions and non-upwelling regions ( Figures 2A, B ), the temperature difference between S S T ¯ r e g E   ( S S T ¯ r e g W ) and S S T ¯ e a s t   ( S S T ¯ w e s t ) can be used to indicate the change of upwelling strength to some extent. A big temperature difference indicates a strong upwelling. As shown in Figure 2C , the increasing SST difference suggests that both the Qiongdong and Qiongxi upwelling are likely to be strong from 2003 to 2021, especially the period after 2013 when the average SST difference is larger than 1.2 °C for the two upwelling regions.

In addition, the PC values of EOF3 mode have a slight declining trend and the values are almost negative after 2013 ( Figure 5C ), which means that the ES upwelling enhances and brings more cold water. Xie et al. (2016) found a weakening Qiongdong upwelling using the index of SST temperature between non-upwelling region and upwelling region during 1982-2012. Our results also show a decreasing SST difference and weakening Qiongdong upwelling before 2013 ( Figure 2C ), which is consistent with their results. However, we found an enhanced trend of the upwelling after 2013 ( Figure 2C ). Although the upwelling becomes strong, the SST in the upwelling regions become warm and the areas of low-temperature water decrease ( Figure 2 ). It meant that the enhanced upwelling does not seem to compete with the overall SST warming in the northwest area of the South China Sea ( Figure 5A ). The overall SST warming in the northwest area of the South China Sea is responsible to the warming trend of SST in upwelling regions around Hainan Island. However, enhanced upwelling slows the warming trend of SST in the upwelling regions compared with other sea regions ( Figures 2A, B ), which is beneficial to the marine ecosystem and fisheries activities.

In addition to the warming trend, SST in upwelling regions also showed evident interannual variations ( Figure 4B and Figures 5B, C ). Wind is an important factor in inducing coastal upwelling. Previous studies have shown that Ekman transport caused by the southwest monsoon wind and Ekman suction caused by the local wind stress curl are two dynamic mechanisms for the generation of Qiongdong upwelling (Jing et al., 2009). For the Qiongxi upwelling, the convergence caused by the southwesterly wind hinders the formation of upwelling (Wang et al., 2015b). To explore the impact of wind on the interannual variation of upwelling, we calculated the wind stress and wind stress curl in the ES (110–112°E, 18–19.64°N), NES (110–112°E, 19.64–21.50°N), and WS (108–109°E, 18–19.64°N) upwelling regions, as well as the along-coast wind stress in the ES upwelling region (Section 2.1). Correlation analysis showed that the interannual variations in wind stress and along-coast wind stress do not correlate with the dSST anomaly, which indicates that the wind stress itself may not be the main factor affecting the interannual variation in upwelling, which is consistent with the conclusion of Xie et al. (2016).

The time series of wind stress curl in the three upwelling regions is shown in Figure 6 . The wind stress curl is positive in the ES and NES upwelling regions in summer, whereas it is negative in the WS upwelling region (except in 2012 and 2018). In terms of Ekman suction, the positive wind stress curl over the ES and NES upwelling regions is beneficial to the generation of upwelling, whereas the negative wind stress curl over the WS upwelling region can inhibit upwelling. Correlation analysis showed that the wind stress curl in the ES upwelling region and PC3 had a significant negative correlation of -0.63 (p< 0.01), which indicates that the wind stress curl can be an important factor for interannual changes in SST over the ES upwelling region and part of the NES upwelling region. Meanwhile, the wind stress curl had a significant negative correlation with the interannual variation in dSST anomalies (r = -0.54, p< 0.05) in the ES upwelling region, indicating that the interannual variation in ES upwelling is closely affected by the wind stress curl. However, the wind stress curls in the NES and WS upwelling regions did not correlate with the interannual SST variation in the corresponding upwelling regions.

The NES upwelling is caused by both alongshore wind stresses and large-scale circulation in the South China Sea. The large-scale circulation flowing northeastward along the east coast transports cold water from the south and turns offshore off the northeast coast, uplifting the water mass from the bottom layer near the shore (Li et al., 2012; Lin et al., 2016a; Lin et al., 2016b). Lin et al. (2016b) suggested that the uplift of cold bottom water can be limited by strong stratification, which is associated with the prevailing southeasterly wind and fresh water, hindering the occurrence of NES upwelling on the sea surface. Figure 7A shows the relationship between the SST in the ES and NES upwelling regions. In most years, the interannual SST variations in the two upwelling regions were consistent (near the 1:1 line). However, the SST of the NES upwelling region (29.00°C) was considerably higher than that of the ES upwelling region (27.89°C) in 2018, suggesting the NES upwelling intensity was markedly weak. We further analyzed the interannual variation in total precipitation and wind speed over the NES upwelling region ( Figure 7B ). In 2018, the precipitation in summer is the highest (0.57 mm hr^-1) and weak southeasterly wind (1.82 m s^-1) prevails. Both are conducive to the southward expansion of fresh water to the NES upwelling region, thus forming a strong density stratification that hinders the uplift of bottom cold water, causing a weak NES upwelling and a high SST in the region.

The formation of WS upwelling may be related to the summertime bottom cold water mass in the Gulf of Tonkin (Gao et al., 2014). In summer, there is a strong horizontal pressure gradient force between the stratification zone in the center Gulf of Tonkin and the vertical mixing zone near the coast of Hainan Island, which can generate upwelling across the tidal front from the stratification region to the mixing zone. Upwelling occurs at the edge of the bottom cold water mass and brings the bottom cold water below the thermocline to the sea surface, which is similar to edge upwelling around the Yellow Sea bottom cold mass in summer (Zhu et al., 2018). Previous studies have suggested that the bottom cold water mass in the central Gulf of Tonkin during summer is closely related to the retention of local cold water during winter (Gao et al., 2017). The sea region of 18–20.5°N and 106.5–108.5°E is selected to represent the central Gulf of Tonkin. Since the water volume is well-mixed in winter, we calculated the correlation coefficient between the interannual variation of summer SST in WS upwelling region and the winter SST in the central Gulf of Tonkin, and found that they were significantly positively correlated (r = 0.53, p< 0.05) ( Figure 8 ). This indicates that the interannual variation in the WS upwelling is associated with the water temperature in the central Gulf of Tonkin. As for the Qiongdong upwelling, it should be noted that the relevant indicators in the ES and NES upwelling regions do not correlate with the interannual variation in water temperature in winter.

The previous discussion suggested that SST in winter is an important factor affecting the interannual variation in summer WS upwelling. We further examined the impact of wintertime atmospheric conditions. The winter SST in the central Gulf of Tonkin was associated with the local air temperature (r = 0.87, p< 0.01). The correlation analysis showed that the air-sea heat flux in the central Gulf of Tonkin is correlated with interannual variations in local air temperature (r = 0.68, p< 0.01) and local SST (r = 0.46, p = 0.06) in winter ( Figure 8 ). This indicates that the winter air temperature associated with air-sea heat flux is likely to affect the SST in the central Gulf of Tonkin in winter, and then affect the interannual variation of WS upwelling in summer.