Seasonal Shift of Storm Surges in the Yangtze Estuary, China

Abstract

This study explores storm floods in the Yangtze Estuary to investigate how extreme sea levels and storm surges change in the context of global warming. Previous studies focused on the long-term variations in amplitude or frequency of storm surges, with limited research conducted on the timing of extreme storm surge events. Based on the methods of non-stationary extreme value theory, we explored the last 33-year tidal levels at Xuliujing Station and found that the annual extreme water level has exhibited a slight downward trend, which is directly attributed to the decrease in mean sea level resulting from reduced upstream river flow. The storm surge season of the Yangtze Estuary experienced a significant lag in the period 2005–2018, which is not restricted to the Yangtze Estuary but is rather a large-scale climate characteristic of a broad oceanic region. The reason for this shift is the sustained increase in the intensity of the western Pacific subtropical high in the last 15 years, leading to the appearance of low-pressure channels in the East China Sea in September and October and thus causing more typhoons to enter the East China Sea during the later period of the storm surge season.

1. Introduction

In the context of climate change, extreme weather and climate events occur frequently [1,2]. Since the 1970s, there has been a significant increase in the proportion of super typhoons, as well as the number of strong typhoons and associated rainfall making landfall in the coastal areas of East Asia and Southeast Asia, attributed to the warming global ocean [3,4,5]. According to the Sixth Comprehensive Assessment Report of the IPCC, global surface temperatures have risen by 1.1 °C from the period of 1850–1890 to the period of 2011–2020 [6]. With each degree Celsius of temperature increase, there is a potential increase in the proportion of super typhoons by 25% to 30% [7]. By 2100, there may be a 14% increase in the average typhoon intensity in the Northwest Pacific [8]. Estuaries are susceptible to a sudden rise in water levels induced by intense typhoons and abnormal heavy rains, significantly elevating the risk of flooding in low-lying areas [9], thereby posing a major threat to coastal society and urban development [10]. Therefore, it is crucial to investigate the trend of changes in extreme sea levels and storm surges in the warming climate, which has become a prominent topic in current research on coastal disaster prevention and mitigation in estuarine areas.

Sea level rise (SLR) is an enduring characteristic of climate change [11]. SLR has been the primary driver of extreme sea level events since the 1960s, based on available tide gauge records [12]. By 2100, the global mean SLR is projected to reach 0.61–1.10 m under the RCP8.5 scenario [13]. SLR will elevate estuarine base water levels, potentially causing beach erosion and shoreline recession, thereby jeopardizing basic coastal hydraulic facilities [14,15]. Under the influence of storm surges, the rising sea level has contributed to an increased frequency of extreme tidal levels worldwide, exacerbating the risk of coastal inundation [16]. However, storm surge heights have not exhibited a consistent trend on a global scale. Mawdsley and Haigh [17] found that significant trends in storm surge were observed at only 13% of the sites among 220 global tidal stations, with approximately an equal number exhibiting positive and negative trends. Despite storm surges contributing significantly to extreme tidal events [12], research conducted in the North Atlantic region has found no clear trends in surge changes, whether in Europe or along the North American coast [11]. However, in the South China Sea region, there has been a noticeable increase in the frequency, duration, and intensity of storm surges since the late 1980s [18]. This suggests that global trends in storm surge variations are inconsistent and exhibit a distinct regional feature.

In a warming climate, tropical cyclones also show distinct long-period characteristics. Daloz and Camargo [19] found that tropical cyclones around the world have shown a trend of migrating poleward in the past few decades. Niu [20] suggested a significant increase in the number and intensity of tropical cyclones in the Western Pacific, along with their northward expansion in terms of impact coverage. Under the RCP8.5 scenario, both the genesis locations and the locations of maximum intensity of typhoons are projected to shift northward [21,22]. The activity range of tropical cyclones has expanded from tropical to subtropical regions. Guo and Tan [23] indicated that the eastward retreat and northward shift of the western Pacific subtropical high, combined with La Niña-like sea surface temperature conditions, constitute the main drivers behind the northward shift of typhoon tracks in the Western Pacific. Additionally, typhoon activity increased significantly in September and October in East Asia from 1981 to 2019 [24]. Liang [25] found that the seasonal changes in super typhoon intensity lag behind the scale seasonal changes by two months.

The escalation in extreme sea level events is primarily attributed to two factors: sea level rise and changes in storm surge events. Climate change has led to regional variations in tropical cyclone patterns globally, while the trend in storm surge events is not consistently apparent. Moreover, prior studies solely concentrated on the long-term variations of surge amplitude or frequency; investigations into changes in the timing of extreme surge events have rarely been conducted. Reinert et al. [26] discovered that the storm surge season along the European Atlantic coast is manifesting an earlier onset, yet the cause of this phenomenon remains unclear. Roustan et al. [27] further found that the shift of the storm surge season along the Northeast Atlantic coasts is related to the North Atlantic Oscillation. In the context of climate change, the characteristics of long-term variation in surge heights and surge timing are rarely studied. The primary objective of this paper is to investigate the long-term trend of storm surge heights and surge season, along with the underlying mechanisms, utilizing long-term tidal gauge records and reanalyzing meteorological data, incorporating methods from extreme value theory.

2. Data and Methods

2.1. Data Description

2.1.1. Water Level and River Flow Data

To investigate the historical long-term changes in storm surges in the Yangtze Estuary, this study utilizes river flow data from the Datong Hydrological Station and tidal gauge records from three stations near the Yangtze Estuary: Xuliujing, Yanglin, and Gongqingwei (

Table 1. Yangtze Estuary and Datong hydrographic station data information table.

Station Name	Lon (°)	Lat (°)	Measurement Item (s)	Period of Measurement
Xuliujing	120.93	31.75	Hourly tide level	1988–2020
Yanglin	121.27	31.58	High/Low tide level	1988–2020
			Hourly tide level	2013–2020
Gongqingwei	121.84	31.38	High/Low tide level	1988–2020
			High/Low tide level	2017–2020
Datong	117.68	30.78	Daily average flow rate	1956–2020
			Hourly flow rate	2003–2020

and Figure 1). Datong station in Anhui province, China, is situated approximately 624 km from the river mouth, representing the tidal limit of the Yangtze River. The Datong flow rate data primarily comprise two components: the daily average discharge and the hourly discharge. The daily average discharge is extracted from the hydrological yearbook, covering the period from 1954 to 2003, while the hourly discharge represents the measured flow and spans from 2003 to 2020. The tidal level data at the mouth of the Yangtze River comprise two components: hourly water levels and daily high/low tide data. The hourly water level data at Xuliujing Station spans from 1988 to 2020, while the data at Yanglin and Gongqingwei Stations span from 2013 to 2020 and from 2017 to 2020, respectively. The tidal datum is referenced to the Wusong Datum. Before utilization, the tidal data undergo cleaning, validation, and interpolation processes to address any data gaps owing to the manual transcription of early tidal data.

Among the three stations, Xuliujing Station serves as the starting point for the lower section of the Yangtze Estuary. It exhibits a slightly higher sensitivity to the impact of upstream runoff and constitutes the primary focus of this study. Gongqingwei Station is situated on the southeastern side of Hengsha Island, adjacent to Chongming Island, and faces the open sea. The influence of upstream runoff on the tide level is minimal, amounting to just a few centimeters, and thus can be disregarded. Yanglin Station is positioned between the other two stations. Considering Gongqingwei Station as the reference point, the distances from Yanglin and Xuliujing Stations along the river axis are 55.7 km and 96.8 km, respectively (Figure 1).

2.1.2. Western Pacific Subtropical High Data

The western Pacific subtropical high (WPSH) is one of the primary circulation systems impacting weather and climate in East Asia. Its extent is represented on the 500 hPa upper-air chart by the 588 dagpm contour encircling the region from 110° E to 180° E. Two indices of WPSH are applied in this paper: the intensity index and the western boundary index. The WPSH intensity index is defined as the sum of the product of the area and the difference between the geopotential height value of the region and 587 dagpm, representing the “volume” of the WPSH over 588 dagpm. The region is defined as the area above 10° N and within the range of 110° E to 180° E on the 500 hPa geopotential height field, where all geopotential heights are greater than or equal to 588 dagpm. The western boundary index is the minimal longitude of the 588 diagram contour in the region between 90° E and 180° E, indicating that its minimum value is 90° E [28]. The monthly data for the WPSH intensity index and western boundary index in this paper, covering the period from 1951 to 2020, were obtained from the National Climate Center of the China Meteorological Administration. (http://cmdp.ncc-cma.net/Monitoring/cn_stp_wpshp.php?wpsh_elem=wpsh_GD, accessed on 2 February 2024). The grid data for WPSH were obtained from NCEP/NCAR, with a spatial resolution of 2.5° × 2.5° (https://www.psl.noaa.gov/data/gridded/data.ncep.reanalysis.html, accessed on 2 February 2024).

2.2. Methods for Studying Non-Stationary Time Series

When applying the classical extreme value theory to investigate problems, it is typically assumed that the variables ($X_1,X_2,\dots X_n$) under study follow an independent and identically distributed (IID) distribution. That is, statistical parameters such as the mean and variance of the distribution remain constant across all variables and are not subject to trends, cycles, random walks, or their combination. The condition of “independence” is generally not easy to satisfy. In many natural phenomena, current values ($X_i$) often depend on recent values ($X_{i-1},X_{i-2}\dots$) from the past. For example, there is always a certain correlation between rainfall on two consecutive days. Apart from randomness, many natural phenomena show systematic properties such as trend, seasonality, or periodicity. Such time series are non-stationary time sequences. The parameters of the non-stationary time sequence model are not fixed but change with time. The distribution function of a non-stationary time series is similar to the data-independent distribution function, but its parameters are functions of the covariate (usually t). According to the study by Coles [29], random variables obey the generalized extreme value distribution (GEV), $x_t~GEV(\mu (t),\sigma (t),\xi (t))$, and its cumulative distribution function is:

$$G(x_t;\mu (t),\sigma (t),\xi (t))=\exp \left\{-{\left(1+\xi (t)\frac{x_t-\mu (t)}{\sigma (t)}\right)}^{-1/\xi (t)}\right\}$$

(1)

where $\mu (t)$ is the location parameter, $\sigma (t)>0$ is the scale parameter, and $\xi (t)$ is the shape parameter. Generally speaking, the shape parameter is difficult to estimate and exhibits limited change over time [30,31], which is always assumed to be constant. That is, $\xi (t)=\xi$, and its value is determined at the first estimate. In particular, the above distribution is widely known as the Gumbel distribution when $\xi (t)=0$ [29].

The location parameter $\mu$ in the GEV distribution of non-stationary time series is defined as:

$$\mu (t)={\mu }_0+{\mu }_{1}t+{\mu }_2\cos (\omega t)+{\mu }_3\sin (\omega t)+{\mu }_4\mathrm{SOI}(t)$$

(2)

where ${\mu }_1$ represents the linear trend (e.g., the long-term trend of storm surge); the harmonic oscillation terms (the coefficient ${\mu }_2,{\mu }_3$ and $\omega =2\pi /yr$) represent the seasonal characteristics of storm surge [32]; the Southern Oscillation Index (SOI) is introduced to examine the effect of climate on extreme storm surges, with the coefficient ${\mu }_4$. The SOI data are from the National Oceanic and Atmospheric Administration (https://www.cpc.ncep.noaa.gov/data/indices/soi, accessed on 2 February 2024).

The scale parameter $\sigma$ ($\sigma >0$) is defined in the form of the location parameter $\mu$:

$$\sigma (t)={\sigma }_0+{\sigma }_{1}t+{\sigma }_2\cos (\omega t)+{\sigma }_3\sin (\omega t)+{\sigma }_4\mathrm{SOI}(t)$$

(3)

The GEV model was determined using the maximum likelihood method to obtain the parameters, and the standard error of each parameter estimate was calculated via the delta method described by Coles [29]. It should be noted that the GEV model requires the maximum height and the exact time of the storm surge as input. At the same time, the skew surge is a variable calculated based on the entire tidal cycle and does not have any specific time-point attribute. Therefore, the time of the skew surge is defined as the time of the highest water level in a tidal cycle.

Three methods were employed in this paper to investigate the changes in the seasonality of extreme surge levels. The first method is “sliding window analysis” (method 1): We applied the non-stationary GEV model to every 20-year window and estimated the model parameters. In each window, the phase of the seasonal oscillation is analyzed, which means the timing of the maximum of ${\mu }_2\cos (\omega t)+{\mu }_3\sin (\omega t)$, cf. Equation (2). This represents the timing of extreme surge events in the year. The method of “sliding window analysis” requires at least 50% data completeness. That is, every continuous period of 20 years should contain 120 monthly maxima or more; otherwise, the estimation will terminate [26].

The second method is “inverse distance weighted”: the timing of extreme surge events in a year is calculated using the timing and height of storm surges directly and is defined as:

$$e={\displaystyle \sum _{i=1}^n{t}_i}{S}_i/{\displaystyle \sum _{i=1}^n{S}_i}$$

(4)

where $n$ is the number of storm surges that occur each year; $S_i$ and $t_i$ ($t_i$ is denoted as the time in days of a year) are the height and time of the $i^{th}$ surge, respectively. The method of “inverse distance weighted” assumes that the middle moment of surge season in each year is concentrated in a specific period when violent surges occur. Given that only effective surge height and time are involved in the calculation, there is no need to consider seasonal oscillations.

The third method is “monthly analysis”: The monthly maxima dataset derived from a long time series can be regrouped into 12 new ones by month, to which we fitted the GEV model representatively and obtained 12 sets of parameters. The data in each new time series (monthly maximums) come from the same month in different years, so there is no need for the annual cycle since only a single month is considered in each model.

Therefore, in the third method the parameters are set as ${\mu }_2={\mu }_3=0,{\sigma }_2={\sigma }_3=0$. Similarly, the Southern Oscillation Index in the same month of different years has limited difference numerical terms, so the GEV model estimation parameters of monthly analysis can be simplified as linear:

$${\mu }^{(m)}(t)={\mu }_0^{(m)}+{\mu }_1^{(m)}t,{\sigma }^{(m)}(t)={\sigma }_0^{(m)}+{\sigma }_1^{(m)}t,{\xi }^{(m)},m\in \{\mathrm{Jan},\mathrm{Feb},,..,\mathrm{Dec}\}$$

(5)

The third method has more parameters and can consequently show more changes besides the timing shifts of the extreme surge season [26]. Note that method 3 has a model for each month, so it is more suitable for estimating variables containing the maximum value in each month (e.g., discharge and precipitation). Given that typhoons affecting the middle and lower reaches of the Yangtze River generally occur from June to September, method 3 may not be able to fully leverage its advantages in the study of the seasonal shift of surges at the Yangtze Estuary.

3. Results

3.1. Long-Term Skew Surges Calculation

A storm surge is an abnormal rise or fall of seawater levels caused by atmospheric disturbances and is generally defined in two ways. The first way to define storm surge is the difference between the observed tidal level during a storm and the astronomical tidal level at the same time, which is usually called non-tidal residual in research. The second way is the difference between the storm high tide level and the astronomical high tide level within a tidal cycle (i.e., the ability of meteorological elements to raise the high tide level within the current tidal cycle), which is generally called “skew surge” [33]. Non-tidal residual is used to describe the real-time elevation magnitude of the water level above the astronomical tidal level during a storm, while skew surge is the lifting amplitude of the typhoon to the high water level in a tidal cycle. Coastal defense structures are primarily designed based on the extreme sea levels. Compared to the first form of storm surge (non-tidal residual), skew surge can better depict the risk that typhoons pose to coastal protection. Additionally, the Yangtze Estuary features strong tidal asymmetry, where storm surge calculated in the form of skew surge is easily affected by the phase deviation of the storm tide, and the result is usually larger. Therefore, it is more appropriate to use skew surge as an indicator of storm surge in the Yangtze Estuary from the perspective of disaster prevention and reduction.

To obtain the long-term skew surge from the 33-year observed estuarine water level, the harmonic analysis method is used to hindcast the astronomical tidal level, and the unreasonable result of the skew surge is corrected according to the storm time. According to classical harmonic analysis methods, tides are considered to be the superposition of sinusoidal constituents with different periods. It is generally assumed that the harmonic constants and mean sea level remain constant [34]. However, in the context of climate change and sea-level rise, the harmonic constants are not fixed but change over a long-term period. Additionally, in estuarine areas, the “mean water level” can vary even within a day or a year due to the nonlinear effects of upstream runoff. The variability in the “average water level” results in the non-stationarity of both the amplitude and phase lag, so a method of NS_TIDE [35] for harmonic analysis of non-stationary signals is applied in this paper. NS_TIDE was developed based on T_TIDE [34], taking into account the nonlinear effects of upstream flow on tides. Gongqingwei Station is located south of Hengsha Island, adjacent to the sea, and experiences minimal upstream flow impact, which can be ignored. Therefore, Gongqingwei Station is employed as the offshore reference station for input data in NS_TIDE. Figure 2 shows the comparative predicted tidal level analysis results between NS_TIDE and T_TIDE.

There is a clear positive correlation between the observed tidal elevation at Xuliujing Station and the Datong flow rate. As flow rates increase, the tidal level ensemble at Xuliujing Station rises, while a decrease in flow rates results in a corresponding decline in tidal level. Examining the root mean square errors, we find that NS_TIDE has an error of merely 0.18 m, while T_TIDE’s error stands at 0.19 m, indicating a marginal disparity between the two methodologies. Notably, despite occasional small errors in some periods, both techniques excel in predicting the astronomical tidal levels at Xuliujing Station, as illustrated in Figure 2a. However, T_TIDE assumes that the sea level is constant and aims to minimize the error using the least squares method for harmonic analysis. This inherent characteristic leads to T_TIDE’s harmonic water level results converging towards the “mean,” resulting in underestimation during flood periods and overestimation during dry periods. The comparison of water level predictions during flood periods between the two methods is presented in Figure 2b. When the flow rate is relatively high, approximately 70,000 m³/s, T_TIDE’s calculations display a slight underestimation of 0.30–0.40 m (depicted by the green solid line), while NS_TIDE’s calculations exhibit only a modest overestimation of around 0.10 m (illustrated by the red solid line), which is more consistent with the observed tidal levels. The comparison of water level predictions during dry periods between the two methods is presented in Figure 2c. When the flow rate is relatively low, approximately 15,000 m³/s, Xuliujing Station experiences overall low water levels. T_TIDE’s calculations indicate a modest overestimation of around 0.20 m (depicted by the green solid line), while NS_TIDE’s errors are less than 0.10 m (illustrated by the red solid line). Evidently, in scenarios where upstream flow rates are either notably high or low compared to the annual average, the tidal level prediction results of NS_TIDE are superior.

The primary objective of this section is to separate astronomical tides from observed tides and calculate storm surges. It focuses more on the effectiveness of harmonic analysis during typhoon events than demanding high overall accuracy for annual harmonic analysis. Strong typhoons that impact the Yangtze Estuary typically occur during the summer and autumn, coinciding with the Pacific warm and humid airflow brought by the summer southeast monsoon. This airflow often results in prolonged and extensive rainfall in the middle and lower reaches of the Yangtze River, thus leading to riverine flooding. NS_TIDE provides more accurate results during flood periods. Therefore, compared to T_TIDE, using NS_TIDE for astronomical tide calculation and storm surge estimation at Xuliujing Station yields higher accuracy and more reliable results. Furthermore, the encounter of storm surges and flooding represents a more extreme and potentially more devastating scenario, highlighting the necessity of employing non-stationary tidal harmonic analysis for astronomical tide prediction.

It should be noted that in certain scenarios, the harmonic analysis may introduce minor phase deviations in the reported astronomical tides due to the incomplete separation of subtidal tidal constituents and other high-frequency tidal constituents [11]. The phase discrepancy in estuarine tidal constituents significantly affects the computation of storm surge (non-tidal residuals). In contrast, employing skew surges as an indicator for estuarine storm surges effectively mitigates the impact of phase shift-induced errors.

A time series of astronomical tides spanning 33 years from 1988 to 2020 can be extracted using harmonic analysis with NS_TIDE. By calculating the difference between the highest astronomical tide and the corresponding observed highest tidal level within each tidal cycle, a time series of skew surges can be obtained. A more conservative approach, referred to as the ISYE-R500 method [9], is employed for the identification of storm tracks and verification of skew surges. Figure 3 presents the time comparison between typhoons impacting the Yangtze Estuary and skew surges at the Xuliujiang Station. The gray boxes in the figure indicate the occurrence times of typhoons throughout the year, while the solid dots represent the times of skew surges at the Xuliujiang Station, a larger dot with darker colors indicating larger values. Notably, the timing of surges at the Yangtze Estuary aligns closely with the times when typhoons impact the region. Each circle is encompassed by a gray box, although there are a few instances where circles are not enclosed. This indicates that ISYE-R500 can identify all typhoons that generate surges at the Yangtze Estuary, but there is a probability of incorrectly categorizing typhoons that do not produce surges, with an accuracy rate of approximately 81%. Consequently, in the absence of observed tidal-level references, using only typhoon track data and the ISYE-R500 method to identify typhoons leading to storm surges in the Yangtze Estuary is a reasonably conservative approach. From the figure, it is evident that storm surges in the Yangtze Estuary predominantly occur between August 1st and October 1st each year, accounting for a significant 67% of the total occurrences. Over the past 33 years, the earliest peak storm surge occurred on May 18th (in 2006, Typhoon Pearl) with a surge value of 0.25 m, while the latest was observed on November 24th (in 2019, Typhoon Phoenix) with a surge value of 0.60 m. The largest peak storm surge, measuring 1.48 m, was recorded in 1997 (during Typhoon Winnie).

3.2. Trend Component of Extreme Tidal Levels

At first, we investigated long-term trends in our tide gauge record to understand how the data generally evolved. We calculated the annual MSL, the annual 99th percentile of the full sea level (including MSL, tide, and surge), and the annual 99th percentile of surge levels. To investigate how estuarine water levels evolve in the context of climate change, an analysis of the long-term trends in the annual MSL and extreme tidal levels at the estuary is conducted. Figure 4 displays the annual variability in the 99th percentile tidal levels at the Xuliujing Station from 1988 to 2020. Notably, there is significant variability in annual extreme tidal levels, with a standard deviation of 0.10 m. The lowest recorded extreme tidal level at the Xuliujing Station was 4.51 m (in 2006), while the highest extreme tidal level was 4.97 m (in 1998). When considering the long-term trends, the Xuliujing Station exhibits a mild decline in annual extreme tidal levels, with a rate of −1.60 mm/a. This decline in extreme tidal levels may be associated with changes in MSL.

Figure 5 illustrates the long-term variations in MSL at the Xuliujing and Yanglin Stations. Evidently, the interannual changes in MSL at the Xuliujing Station are significant. Excluding the peak value of 2.93 m recorded in 1998 and the minimum value of 2.62 m in 2011, the MSL predominantly fluctuates between 2.70 and 2.85 m. Concerning long-term trends, the annual MSL at the Xuliujing Station undergoes a gradual decline with a rate of −0.7 mm/a. This suggests that the descending MSL at Xuliujing Station is one of the contributing factors to the decline in extreme tidal levels.

Research on the long-term patterns of extreme tidal levels has revealed that, globally, most tide gauge stations show an increase in extreme tidal levels. However, this trend diminishes significantly when adjusted for changes in mean sea level [11,32]. Hence, changes in extreme tidal levels are believed to be primarily attributed to the rise in mean sea levels, aligning with the findings of this study. Findings from Sweet and Park’s research [14] additionally suggest that long-term sea level rise is a significant driving factor for increased flood hazards along the U.S. coastline. Despite the slight decline in mean sea level at the Xuliujing Station, this trend contradicts the global pattern of rising temperatures, melting ice caps, and a general sea level increase. As an example, the Dajishan tide gauge station, located outside the Yangtze Estuary, has observed an average sea level rise rate of 3.15 mm/a over the last 42 years [36].

The deceleration in the mean sea level rise rate within the Yangtze Estuary and its counter-trend to the sea level rise in the open sea can be attributed to changes in upstream discharge. Illustrated in Figure 5b, the long-term trend of mean sea level at the Yanglin Station (situated downstream of the Xuliujing Station) is also upward, at a rate of 2.34 mm/a, which is approximately 30% lower than the rate of sea level rise in the open sea. Figure 6 shows the relationship between the annual mean water level at the Xuliujing Station and the annual mean discharge at the Datong Station.

Analysis of this graph reveals that in the past 30 years, the annual mean discharge at the Datong station has consistently decreased, with an annual decrease rate of −460 m³/s, aligning with the trend in the annual mean water level at the Xuliujing Station. Interannual variations in MSL and Datong discharge trends exhibit a close correlation. During years of increased Datong river flow, the annual MSL at Xuliujing correspondingly rises compared to the previous year, whereas decreased Datong river flow leads to a decline in the annual mean water level at Xuliujing, emphasizing a strong correlation between the two ($r=0.823$).

3.3. Trend Component of Extreme Skew Surges

Annual maximum skew surges were extracted at the Xuliujing Station, revealing that, except for the extreme skew surge of 1.48 m during Typhoon “Winnie” in 1997, over 90% of these surges are below 1.0 m (Figure 7). From a long-term perspective, the annual maximum of peak storm surges exhibits a declining trend, with an average decrease rate of 1.65 mm/a. The impact of MSL changes on skew surges is canceled by employing the non-stationary tidal harmonic analysis method for calculating astronomical tides. Therefore, the decrease in annual maximum skew surges is not the result of the declining MSL at Xuliujing Station; climate change is likely the primary contributing factor. We further applied the non-stationary GEV model (Equation (2)) to the “monthly maximum” skew surges from 1988 to 2020 and the estimated location parameter ${\widehat{\mu }}_1=0.1$ mm/a. Although the linear trend is statistically significant ($\alpha =0.05$), its magnitude is very small compared to changes in sea level and extreme tidal levels. Therefore, there has been no significant trend in skew surge at the Xuliujing Station over the past 33 years.

As a secondary consequence of typhoons, the dynamics of storm surges are additionally impacted by large-scale climate features. In the non-stationary GEV model for “monthly maximum” skew surges over the past 33 years, incorporating the Southern Oscillation Index (SOI) as an input reveals that the impact of SOI on estuarine skew surges is statistically significant, but the value is relatively small. A one-unit change in SOI can lead to an approximate 0.3 cm variation in the location parameter (${\widehat{\mu }}_4=0.30\pm 0.21$). Furthermore, the findings ${\widehat{\mu }}_4>0$ indicate that the probability of typhoon-induced surges is greater in positive SOI years, aligning with the previous studies [37,38]. In El Niño years, where typhoon activity decreases (with a negative SOI), and La Niña years, where typhoon activity increases (with a positive SOI), the corresponding likelihood of triggering extreme surges in the Yangtze Estuary also increases. Considering the relatively modest impact of SOI on storm surges, we conducted a re-estimation (let ${\mu }_4={\sigma }_4=0$) and demonstrated that including or excluding the influence of SOI has minimal impact on the examination of seasonal variations in surges. Therefore, in the subsequent investigation of seasonality, the influence of SOI is disregarded.

3.4. Seasonal Characteristics of Skew Surges

Storm surges induced by typhoons in the Yangtze Estuary exhibit pronounced seasonal characteristics. As typhoons in the western Pacific commonly arise during the summer and autumn, the influence of typhoons and their related surges in the Yangtze Estuary is predominantly observed between June and October annually, constituting over 90.5% of occurrences. No evidence of surges was observed between January and April. In the non-stationary GEV model, the amplitudes representing the interannual oscillation for the location and scale parameters are ${\left({\widehat{\mu }}_2^2+{\widehat{\mu }}_3^2\right)}^{1/2}=1\pm 0.2$ cm and ${\left({\widehat{\sigma }}_2^2+{\widehat{\sigma }}_3^2\right)}^{1/2}=0.7\pm 0.2$ cm, respectively. The phase angles are $\arctan ({\widehat{\mu }}_3/{\widehat{\mu }}_2)=221\pm 11°$ and $\arctan ({\widehat{\sigma }}_3/{\widehat{\sigma }}_2)=196\pm 14°$, respectively. The estimated location parameter $\widehat{\mu }(t)$ reaches its maximum in early August, while the estimated scale parameter $\widehat{\sigma }(t)$ peaks in mid-to-late July, indicating that the two parameters are asynchronous in time (although to a small degree).

Besides the height of the extreme storm surge, changes in the timing of the surge can also be captured in the non-stationary GEV model with a sliding window of a 20-year scale (Method One). Given the omission of the influence of climate indices, the day in a year corresponding to the maximum of the interannual oscillatory component of the location parameter $\mu$ is identified as the central moment of the storm surge season, representing the timing of extreme surges. The timing of extreme surges typically falls between August and September (solid blue line in Figure 8), with a 95% confidence interval covering approximately 2 to 3 weeks. There is an evident trend that, before 2005, the extreme storm surge season consistently shifted earlier, reaching the end of July, while after 2005, there was a noticeable delayed trend in the extreme storm surge season, extending into mid- to late-August by 2018. The orange solid line in Figure 8 illustrates the outcomes of “inverse distance weighted” (Method Two), which underwent an 8th-order Butterworth low-pass filter to capture the intergenerational variations in storm surge seasons. Despite Method Two yielding results approximately 13 days later than Method One on the whole, the trends and amplitudes of the two methods in the shift of storm surge seasons are highly consistent. That is, before 2005, the extreme storm surge season displayed an advancing trend, whereas from 2005 to 2018, it exhibited a discernible lagging trend.

The primary determinants of the overall numerical disparities in calculating the timing of storm surges between the two methods are inherent characteristics. The input data of the non-stationary GEV model with a sliding window are the monthly maximum surge values, resulting in a continuous distribution of extreme storm surges. This distribution is less sensitive to individual large events and is more representative. The maximum surge in typhoon months is caused by typhoons, whereas the maximum surge in the other months is influenced by other meteorological factors and errors in the harmonic analysis. Even though the values are generally less than 5 cm, their inclusion in the GEV distribution still results in a slight shift towards the first half of the year. The “inverse distance weighted” method takes into account only the actual occurrence times and magnitudes of storm surges in a year. It can effectively handle scenarios in which multiple strong typhoon surges occur in a single month, rather than just the monthly maximum. As a result, its outcomes are more closely aligned with the actual occurrence of typhoons throughout the year (July to October), resulting in a shift towards the second half of the year. However, results obtained through this discrete data calculation are more influenced by individual strong typhoon surge events, making them less representative than Method One.

Both methods have their merits, and their results collectively indicate a delayed trend in the extreme storm surge season in the Yangtze Estuary from 2005 to 2018. Therefore, the average of the results from both methods is used to represent the trendline (red dashed line in Figure 8). The combined result reveals that, during 2005–2018, the extreme storm surge season in the Yangtze Estuary was delayed from 7 to 24 August, representing a delay of 17 days.

3.5. Seasonal Characteristics of Skew Surges in Large Areas

Typhoons affecting the Yangtze Estuary generally originate over the expansive ocean area east of the Philippines, traverse the Taiwan–Ryukyu Islands region, and then enter the sea areas of the East China Sea. Influenced by subtropical high-pressure systems, the northward trajectory of these typhoons typically follows a concave path from the southeast to the northwest. To investigate whether the surge season shift at Xuliujing was limited to the Yangtze Estuary or extended to a larger regional sea-air feature during the period from 2005 to 2018, we conducted a harmonic analysis and trend study on hourly sea level data from six tidal stations along the Hong Kong–Kaohsiung–Ryukyu Islands coast. The data were obtained from the Global Extreme Sea Level Analysis (https://gesla787883612.wordpress.com/downloads/, accessed on 2 February 2024). The locations of the tidal stations are shown in Figure 9, and the period of sea level data and the root mean square error (RMSE) of harmonic analysis are shown in

Table 2. The tidal-level data information and harmonic analysis error over a wide area of ocean.

No.	Name	Lon.	Lat.	Period	RMSE (m)
h329	HongKong	114.217	22.300	1986–2020	0.1100
h340	Kaohsiung	120.288	22.615	1980–2016	0.0618
h365	Ishigaki	124.150	24.333	1969–2020	0.0600
h355	Naha	127.667	26.217	1966–2020	0.0598
gs21	Okinawa	127.825	26.179	1975–2019	0.0636
h359	Naze	129.498	28.400	1966–2020	0.0600

. Except for the nearshore station in Hong Kong (Kowloon Bay), which may be affected by nearby topography, resulting in a slightly larger error (RMSE = 0.11 m), the tidal stations near the islands all yield satisfactory results in the harmonic analysis, with errors below 0.065 m.

The surge season of the Yangtze Estuary typically extends from July to October (cf. Figure 8). Using the “monthly analysis” method (Method 3) to fit the independent monthly surge time series for July and October during the period from 2005 to 2018, we obtained the trend parameters ${\mu }_1^m$ of these two months, as shown in Figure 10, where the circular points denote the specific data values and the error bars indicate the 95% confidence intervals. It is evident that, excluding Hong Kong and Kaohsiung station, the trend parameters ${\mu }_1^m$ for July are all negative (depicted by blue bars), whereas the trend parameters ${\mu }_1^m$ for October are all positive (depicted by orange bars) during the period from 2005 to 2018. This indicates that surges in July weakened while those in October strengthened, signifying that the surge season for all stations near the Ryukyu Islands is shifting to October. Despite the fact that the tidal stations at Okinawa Island and Naha are close, the specific numerical results differ. This discrepancy arises from the fact that these two tidal stations are located on the north and south sides of Okinawa Island, so even the same typhoon may produce different surge results. The result of the “monthly analysis” GEV model is sensitive to the surge values. Despite numerical disparities, the trends of surge season in both stations are consistent. Therefore, the change in surge seasons in the near-sea area of the Ryukyu Islands and the lagging trend of surge seasons at the Yangtze Estuary are consistent. This indicates that the delay in surge seasons is not confined solely to the Yangtze Estuary but is a climate change process across a broader sea area. The surges at Hong Kong and Kaohsiung Port increased in July and remained relatively constant in October. Both stations exhibit a slight trend of an earlier surge season, which may be related to their respective locations. Hong Kong and Kaohsiung Port are mainly influenced by typhoons entering the South China Sea, while the sea areas near the Yangtze Estuary and the Ryukyu Islands are predominantly affected by typhoons entering the East China Sea. These two regions are influenced by different typhoon systems.

4. Discussion

4.1. Mechanism of Surge Season Delay

The delay in the storm surge season is attributed to changes in the typhoon season in the East China Sea. During summer, the predominant systems in the mid and low latitudes are the westerly trough-ridge and subtropical high-pressure systems, and their varying strengths dictate the fundamental atmospheric circulation, thus establishing the fundamental movement direction of typhoons. Although numerous factors influence the generation and movement paths of typhoons, this paper primarily examines the reasons for the delay in the typhoon season, focusing on changes in typhoon tracks influenced by variations in the western Pacific subtropical high-pressure system and the westward extension point.

Depending on the strength and position of the subtropical high, typhoon tracks affecting China can be mainly classified into three categories. After their generation over the western Pacific Ocean, typhoons tend to move in a west–northwest direction due to the combined effects of the Coriolis force and the eastward winds on the southern side of the subtropical high. When the subtropical high is strong and the westward extension point is farther west than usual, typhoons tend to move along the southern edge of the subtropical high system toward the Indo-China Peninsula, affecting the South China Sea. This type of typhoon track is referred to as a “Westward” typhoon tack (Figure A4a). When the subtropical high weakens and the westward extension point shifts northward and eastward, typhoons are guided by the northwest winds on the southwestern side of the subtropical high, tending to move towards eastern China. This type of typhoon track is referred to as a “Northwestward” typhoon track (Figure A4b). When the subtropical high further weakens, a low-pressure channel forms near the East China Sea. Typhoons move northward through this channel and, under the influence of the westerlies, turn toward Japan and Korea, forming a “Recurve” typhoon track (Figure A4c). When the subtropical high is weak and the westward extension point is farther east than usual, typhoons are not expected to enter the East China Sea (Figure A4d). These different tracks depend on the varying strength and location of the subtropical high-pressure system and its associated atmospheric conditions.

Figure 11 illustrates the interannual variations in the intensity index of the western Pacific subtropical high-pressure system and the longitude of its westward extension point during the typhoon season (July to October) from 1988 to 2020. The intensity of the subtropical high exhibits an increasing tendency (Figure 11a), with an expanding range leading to a continuous westward shift of the extension point (Figure 11b) in the typhoon season. This trend is more pronounced in September and October. Figure 12 shows the position and extent (588 dagpm characteristic line) of the WPSH in each month of the typhoon season. In September, the extent of WHSP (the blue circles) continuously moved westward over time (except for the year 2001, Figure 12d) and finally broke into two parts. In October, the extent of WHSP (the red circles) continuously moved westward and broke into two parts since 2003 (Figure 12e), leading to a low-pressure channel emerging within the maritime regions of the eastern Philippines and the East China Sea. Therefore, typhoons entered the East China Sea through the low-pressure channel. In 2018, the months of September and October witnessed a similar fragmentation of the subtropical high systems, and the low-pressure channel was wider and shifted further westward (as indicated by the pink area in Figure 12i), enabling more typhoons to enter the East China Sea in the later part of the storm season and thus cause more storm surges. This shift in the typhoon season towards September and October led to the lag of storm surge season in the Yangtze Estuary. At the same time, the number of typhoons entering the South China Sea in September and October decreased, resulting in an early occurrence of storm surges in Hong Kong and Kaohsiung ports (Figure 10).

It is important to note that various factors influence the movement of typhoons, including sea surface temperatures, the westerly trough ridge, and more. Moreover, typhoons themselves can also impact the subtropical high systems, albeit with intricate mechanisms, which are not the primary focus of this study. In this section, we have primarily explored the reasons for the delayed typhoon season in the Yangtze Estuary from a macro perspective, considering the intergenerational intensity and spatiotemporal variations of the subtropical high systems.

4.2. Seasonal Shift of Other Hydrological Elements

The seasonal shifts in typhoon-induced storm surges result from changes in large-scale climate conditions. However, it is worth noting that, in addition to storm surges, a similar seasonal shift is observed in flood events. For instance, in the Alpine–Carpathian region, floods have exhibited a delayed transition from summer to autumn flooding [39]. In contrast, in Western Europe along the Atlantic coast, flood events occurred 3 to 8 days earlier during the period from 1960 to 2010 [40]. In northern-central Sweden, spring floods have advanced by approximately one month [41]. Figure 13 illustrates the seasonal shift in major floods at the Yangtze Estuary from 1956 to 2020. This shift demonstrates notable intergenerational characteristics with distinct stages: Stage one (1956–1975) shows a 5-day delay in the flood season; stage two (1975–1992) exhibits a 2-day advancement in the flood season; stage three (1992–2007) displays a 4-day advancement in the flood season; and stage four (2007–2019) indicates a 4-day advancement in the flood season. Significantly, stage four corresponds to the shift observed in storm surge season (2005–2018). Over the past two decades, the flood season at the Yangtze River’s Datong station has been advancing (in July), while the storm surge season at the Yangtze Estuary has been lagging (in August). This reduced the probability of the two peak events coinciding, thereby lowering the overall risk of compounded flood and storm surge disasters at the estuary.

5. Conclusions

In this investigation, to analyze the historical long-term evolution characteristics of storm surges at the Yangtze Estuary, we first introduce an improved variable-speed triangular interpolation method that interpolates daily high and low tidal levels into hourly tidal levels. Then, we extract the skew surges from the hourly tidal levels at Xuliujing using NS_tide over 33 years. Subsequently, several models of extreme value theory are applied to analyze the trend of extreme tidal levels and the seasonal variation in storm surges at the Yangtze Estuary. Finally, the mechanism of the seasonal shift of storm surges is discussed from a meteorological perspective, and the main conclusions are as follows:

Over the past 33 years, the extreme water level at the XuLiujing Station has exhibited a slight decreasing trend (−1.6 mm/a). This trend is primarily attributed to the decreasing annual mean sea level (−0.7 mm/a), which is driven by a reduction in the discharge at the Datong station. The storm surges have not shown a significant trend, and the influence of the large-scale climate feature, the Southern Oscillation Index (SOI), on them is relatively minor.

The storm surge season at the Yangtze Estuary exhibits a distinct seasonal shift. Between 2005 and 2018, the storm surge season shifted from 7 to 24 August, resulting in a 17-day lag. This seasonal shift is not limited to the Yangtze Estuary but represents a climatic characteristic across a broader oceanic region. The primary cause of this shift is the continuous increase in the intensity of the western Pacific subtropical high over the past 15 years, and a low-pressure channel developed over the East China Sea in September and October, leading to more typhoons entering the East China Sea during the later part of the storm surge season.