Does mean sea level trend mask historical storm surge trend: evidence from tropical cyclones affecting Japan since 1980

Variability in storminess, storm surge, and mean sea level (MSL) can substantially alter coastal hazards associated with extreme sea levels (ESLs). However, the detection and attribution of the past changes in tropical cyclone (TC) activity and associated storm surges are hampered by the inhomogeneous TC records. In this study, we investigate spatiotemporal changes in storm surge levels in Japan from 1980 to 2019, a period when observational platforms including tide gauges and storm records are highly consistent. We find statistical evidence supporting the increase in surge annual maxima in several places including the bay area of Tokyo since 1980. This rate of change is comparable to that observed for MSL rise over the same period. These findings cast doubt on the current hypothesis underlying the flood adaptation plan, which assumes that future surge extremes will remain the same and only considers MSL changes. We demonstrate that the changes in ESL in the last 40 years cannot be explained by the rise of MSL alone. Rather, the northeastward shifting of TC landfall location along with intensifying and widening of TCs, might have altered the likelihood of ESL, including surge extremes. The substantial influence of these TC meteorological variables on surge levels combined with the rise of MSL, suggests that current coastal planning practices including critical heights for flood defenses might be inadequate in the future.


Introduction
Storm surge, a rise in sea level due to a storm including tropical cyclone (TC) and extratropical cyclone, acts as a major driver of coastal flooding and has been responsible for as many as 2.6 million deaths worldwide during the past 200 years (Nicholls 2003). Coastal development results in densely populated low-lying cities that are roughly five times (241 people km −2 ) larger than the global mean (47 people km −2 ; (Neumann et al 2015)), making them vulnerable to devastating societal impacts caused by storm surges. On average, 0.8-1.1 million people are exposed to coastal flooding each year globally (Hinkel et al 2014). Several studies including the Sixth Assessment Report of the Intergovernmental Panel on Climate Change have suggested an increase in extreme sea levels (ESL) by 2100 due to storm surges and mean sea level (MSL) rise with high confidence (Oppenheimer et al 2019, Knutson et al 2020, Muis et al 2020, IPCC 2021. ESL refers to the highest water level reached in a given coastal location during specific time periods. It is primarily a function of storm surge and MSL change, and the other factors include tides and waves (Tebaldi et al 2021). Research by Vitousek et al (2017) indicated that even minor changes in ESL (e.g. +10 cm) could significantly increase flooding frequency in many coastal regions between 2030 and 2050. If effective adaptation and mitigation measures are not implemented, the costs of coastal flooding are projected to be substantial for the global economy (Hinkel et al 2014). The effective implementation of these strategies largely depends on a rigorous understanding of the longterm ESL changes (Calafat et al 2022). While considerable attention has been given to the ESL changes resulting from MSL rise in recent decades (Bromirski et al 2003, Woodworth and Blackman 2004, Church et al 2006, Menéndez and Woodworth 2010, Kirezci et al 2020, Taherkhani et al 2020, the quantification of the contribution of changes in storm surge levels to the ESL change remains a grey area that needs to be clarified (Oppenheimer et al 2019).
The analysis of the long-term changes (e.g. 40 years) in the TC-induced storm surge is limited in many regions due to the unavailability of mid to longterm systematic tide gauge observations. Nonetheless, a few recent observational studies have concluded that storm surge levels are increasing in some regions over time and the trends are mainly influenced by MSL rise (Woodworth and Blackman 2004, Church et al 2006, Menéndez and Woodworth 2010. They interpreted their results as storm surge trends are masked by MSL change. Consequently, previous works on future coastal flood risk assessment, which are crucial for adaptation plans, evaluated the future ESL by adding the projected MSL change onto the present climate's storm surge level statistics (Vitousek et al 2017, Kirezci et al 2020, Taherkhani et al 2020. If this assumption of unchanged future storm surge levels turns out to be invalid (either globally or regionally), adaptation plans will be ineffective (Calafat et al 2022). Due to climate change, the characteristics of TCs are likely to change. Specifically, the proportion of intense TCs is expected to increase in the future from present-day values, despite the decrease of the total number of TCs. These variations can change storm surge activity and thus ESL (Knutson et al 2019, 2020, Bloemendaal et al 2022. Therefore, understanding the compound impacts of changing TC characteristics, storm surge levels, and MSL in altering the likelihood of extreme events (on a global, regional, and local scale) is a pressing matter before making the basis of the coastal adaptation plan final (Lal et al 2012).
The contribution of changes in TC characteristics and MSL to storm surge climatology has been studied separately for different time periods and thus the compound association between them is unclear. Several studies indicated that TCs in recent decades are increasingly approaching coastlines with greater intensity, moving at slower speeds, and decaying at a slower rate after landfall, which likely caused by the effects of climate change (Kossin 2018, Li and Chakraborty 2020, Murakami et al 2020, 2022. Such changes in TC meteorological variables have already made significant implications over storm surge intensity. For instance, a study found that the long-term surge risk in Shanghai is more significantly influenced by changes in TC track and landfall location than by intensity . Additionally, Colle et al (2010), Oey and Chou (2016), and Feng et al (2021) separately showed that TC landfall location, intensity, and translation speed have changed over time and thus modulated surge extremes in their study regions. The variability in TC size also has the potential to change surge levels (Irish et al 2008, Islam and Takagi 2020b, Islam et al 2022, however, it is not being studied thoroughly from a historical perspective-probably because the operational service for TC size monitoring had just begun at the end of the twentieth century. While these previous observational studies have used one to one relationship (e.g. TC track vs. surge; intensity vs. surge) to explain the relationship between changes in storminess and observed surge levels, their joint influence (e.g. change in storm surge as a function of change in TC intensity and size) is still unknown. The deficiencies pose an enormous challenge in understating how storm surge levels modulate ESL under the future climate condition.
Here we provide observational evidence that storm surge levels in Japan have increased since 1980 and the rate of change is comparable to that observed for MSL rise. We utilize homogeneous records of MSL, TC storm surge, best track, and meteorological data, which the Japan Meteorological Agency (JMA; see data and methods) maintains. We reveal a consistent northeastward shift in TC landfall location (both eastern and western Japan) and an increased frequency of strong and large sized TCs, which are likely to be the causes of surge level changes. This study suggests that planners may need to reconsider the current hypothesis regarding the coastal flood adaptation plan, in which future surge extremes will remain the same and only MSL changes are considered.

TC selection
In this study, we utilized the JMA best track data archives (JMA 2021b) from 1980 to 2019. This data includes all TCs that originated in the western North Pacific and made landfall in the main islands of Japan. The best track data acquired during the pre-satellite era (i.e. before 1980) contained heterogeneities and significant uncertainties in the data quality (Chan 2019, Moon et al 2019 and were therefore excluded from our analysis. Although the analyzed period was relatively short compared to the period when historical TC record starts officially (since 1951), the period from 1980 to 2019 represented the longest period covered by the JMA with uniform data quality (JMA 2021b). We followed JMA's definition to identify a landfall TC. A landfall TC is identified when the center of a TC reaches the coastline of the mainland (Honshu, Hokkaido, Kyushu, Shikoku) in Japan, JMA recognizes it as a landfall TC (JMA 2021a). We determined the approximate landfall point where a TC track intersects a coastline using vector data provided by the Geospatial Information Authority of Japan (Geospatial Information Authority of Japan 2015). Our analysis focused on TCs with a maximum sustained wind speed (V max ) of no less than 34 kt during the landfall time frame. This wind intensity threshold is necessary as TCs are not always powerful (i.e. tropical depressions) and may not produce noticeable storm surges. Thus, it is more meaningful to focus on tropical storms (34 kt ⩽ V max ⩽ 63 kt) or stronger TCs (i.e. categories 1, 2 in the Saffir-Simpson Hurricane Wind Scale) when managing disaster risk. It should be noted that extratropical cyclones and non-landfalling TCs that pass near the main islands of Japan can also cause storm surges. However, we considered storm surges induced by landfall TCs only. The justification for this is twofold. First, this study primarily focuses on hurricanes/typhoons or tropical storms, collectively known as TCs, rather than extratropical cyclones, which are winter storms. Second, the inclusion of storm surges caused by nonlandfall TCs would require a separate analysis and interpretation of their distinct characteristics. For instance, discussing the meteorological conditions of TCs during landfall would not be possible, as our study specifically aimed to examine the direct influence of meteorological characteristics during landfall on observed storm surges. Based on the aforementioned criteria, 84 TCs were selected for the analysis, which can be regarded as the most economically damaging TCs in Japan.

Tide gauge data
Here, TC-induced peak storm surge (not skew surge) is defined as a residual water level peak that remains after deducting the MSL and predicted astronomical tide from the observed peak storm tide. We considered JMA (JMA 2022) and Japan Coast Guard (Japan Oceanography Data Center 2021) operated tidal stations to retrieve peak storm surge information for each TC. It is important to note that JMA and the Japan Coast Guard publish their tidal observations only after controlling the quality of the initial records. Therefore, the datasets are free from any abrupt changes such as artificial jumps or shifts. The observed tidal records are hourly, and as a result, the peak surge records used in this study may not necessarily represent the exact peak values. Among the many operational stations, data collection was limited to those that met the following five criteria: (a) falling to the right side of a selected TC track and located within the range of radius of 50 kt wind (R 50 ). Radius of 30 kt wind (R 30 ) was utilized when R 50 was not available during TC landfall time frame. While it is possible to use R 30 instead of R 50 in all cases, it will allow us to select several tidal stations for a selected TC track where noticeable storm surges may not be observed (R 30 is generally very large than R 50 that covers several hundreds of kilometers of coasts from TC landfall location). Data base including minor storm surges can affect our statistical analysis. Therefore, we determined R 50 is reasonably a good primary criteria which has been supported by our previous studies, e.g. Islam et al (2022), ); (b) being located on an open coastline or in a bay attached to the main lands (stations on islands were excluded, which constitutes ∼20% of the total tidal stations); (c) having available JMA-predicted astronomical tide data (JMA 2020c); (d) having the elevations of the observation reference plane and the astronomical tide table reference plane available; and (e) having no data gaps occurred when a TC traversed over the station. Based on these criteria, 28 tidal stations were considered in this study. These tidal stations were divided into four regions, namely eastern (R1: 13), western (R2: 7), southern (R3: 6), and northern (R4: 2) Japan (the locations of the tide gauges are shown in figure  S1 in the SI appendix). JMA provides time-series of annual MSL anomaly for each region (same as R1, R2, R3, and R4) (JMA n.d.), enabling us to remove MSL anomaly from respective tide gauge observations and compare surge records with MSL. It needs to be noted that only two TCs (as per TC selection criteria) made landfall in R4 historically and hence, this region was excluded from storm surge analysis in this study. Considering the tidal station selection criteria, peak storm surge information was recorded for more than one station for each TC selected in this study, resulting in a total of 279 surge events during 1980-2019.

TC meteorological data
We analyzed the JMA best track datasets (JMA 2021b) for selected TCs (n = 84) that included TC central positions, intensities (V max ), sizes (R 50 and R 30 ), and translation speeds. Reconstructed typhoon data for Japan (Kubota et al 2021) was used to collect local maximum wind speed recorded at the nearest weather station or lighthouse during the passage of a TC. It needs to be noted that stations used for recording local maximum wind speed are different from tidal stations used for recording surge levels. Thus, we refer to Kubota et al (2021) for full details on the selection of weather stations used for recording local wind speed. Although we analyzed TC meteorological information available from genesis to death, we specifically compared TC meteorological conditions during the landfall time frame with peak surge records. This is because storm surges tend to be amplified during landfall (Islam et al , 2022. However, the TC characteristics during the landfall time frame are not necessarily the most adverse conditions that cause the largest storm surges on the coasts. For example, storm surges would reach the largest value when the TC track is closer to the bay. In this study, we considered the TC landfall time rather than the time closest to each tidal station as a representative condition for causing a peak surge because (a) TC characteristics (i.e. intensity, size, translation speed) after landfall differ from those over the ocean, and thus, the TC information over land is considered less reliable (Huang et al 2021); (b) although the recent TC best track contains more detailed information about when a TC approaches the land, the JMA best track has historically provided TC information at 6 h intervals, making it difficult to identify the time of the closest approach to a tidal station; (c) in the current dataset, surge information was recorded for more than one tidal station for each TC. Hence, a unique characteristic (e.g. landfalling TC intensity) for each storm can provide a simple basis for statistical analysis.

Storm surge hazard potential index
Previous studies (Zhang et al 2000, Grinsted et al 2012 reported that estimating TC potential impact is worthwhile using a storm surge intensity measure. Here, we used the storm surge hazard potential index (SSHPI; Islam et al , 2023 to statistically quantify the contributions of meteorological factors that modulate surge hazards during 1980-2019 in Japan. The SSHPI incorporates meteorological variables that are sensitive to storm surge, including TC intensity, size, and translation speed. Additionally, the SSHPI considers coastal geometry (open coasts and bays) and regional scale bathymetry. These meteorological and geometrical variables are combined into a single measure that represents the expected surge hazard potential along the coast. The TC best track data, particularly during landfall, were used to calculate the SSHPI for each storm. Large surge index values manifest the surge threat posed by a TC at the time of landfall. The bathymetry data for the target region were obtained from the Japan Oceanographic Data Center (Japan Oceanographic Data Center 2020). The effectiveness of the SSHPI for predicting peak surge hazard potential was discussed in . A brief definition is provided in the SI appendix.

Statistical test
The annual surge maxima were assumed to follow a generalized extreme value (GEV) distribution characterized by three parameters: location (µ), scale (σ), and shape (ξ). The parameter µ is assumed to vary linearly with time, while the remaining two parameters are assumed to be time-invariant. The resulting three parameters were estimated from the observations using the Metropolis-Hasting algorithm, which yields uncertainties in estimation through posterior samples of the parameters. Statistical significance for any trend in annual surge maxima provided in this study is based on the 5%-95% credible interval (CI). The return period (p) for a surge event was calculated from the quantile 1-1/p corresponding to the annual surge height of this event. Note that the GEV distribution was only applied to annual surge maxima. For the other variables (i.e. annual MSL and annual mean latitude of peak storm surges), we conducted visual inspection using a quantile-quantile plot and performed Shapiro-Wilk test at the 5% significance level to check if the distribution followed normality. If the data did not significantly deviate from a normal distribution, ordinary linear regression was applied; otherwise, gamma distribution was adopted where parameter µ assumed to vary linearly with time. The 40 year dataset was divided into two sub-periods, 1980-1999 and 2000-2019, hereafter referred to as P1 and P2, respectively. To support the measurement, the non-parametric Kolmogorov-Smirnov (K-S) test at the 5% significance level was used to determine whether samples from P1 and P2 came from the same distribution.

Observational evidence of storm surge change
The change rates of surge annual maxima averaged over eastern (R1), western (R2), and southern (R3) Japan are shown in figure 1(a) (left bars). Hereafter, 5%-95% CI for any surge annual maxima is denoted by square brackets. An overall increase in the surge magnitude of 11.9 mm yr −1 [6.4, 14.1] is evident for 1980-2019. In each region, this change rate varies since storm surge hazards are not spatially homogeneous. While R1 exhibits a maximum change rate of 14.4 mm yr −1 [9.4, 25.3] followed by 9.1 mm yr −1 [0.2, 15.2] in R2, there is no significant change in R3. A more detailed outlook is plotted for the coastal observation stations in figure 1(b). It depicts the substantial difference in the storm surge climatology between eastern and southern Japan. Significant increasing tendencies are noticeable around major coastal cities, including Tokyo (8.47 mm yr −1 [4.1, 11.7]; in R1), Yokohama (7.14 mm yr −1 [4.87, 9.3]; in R1), and Kochi (5.52 mm yr −1 [3.4, 7.72]; in R2). The interquartile ranges associated with the change estimates at each observed station are shown in figure  S2 (see in SI appendix). We further analyze the distribution of storm surge events in Japan between two sub-periods: 1980-1999 (P1) and 2000-2019 (P2) ( figure 1(c)). It reveals a clear shift towards greater annual maxima from P1 to P2 and the two distributions are statistically different at the 5% significance level (K-S test statistic: 0.59; p < 0.05). For example, the occurrence frequency of surge annual maxima greater than 100 cm in P2 is twice as large as that in P1. In this study, each tidal station exhibits different variability in surge levels (e.g. figure 1(b)), nonetheless, combined probability density function shown in figure 1(c) is still useful to demonstrate the overall trend in Japan.
Next, we compare surge change rates (figure 1(a); left bars) with the annual change rates of the MSL during 1980-2019 (figure 1(a); right bars). Here, MSL followed normal distribution (Shapiro-Wilk test statistic: 0.98; p > 0.05), thus ordinary linear regression was applied to estimate change rate (see figure S3 (a) in SI appendix). While the surge to MSL change ratio is more than four times in R1 and R2, the MSL rise in R3 (3.5 mm yr −1 [3.3, 3.6]) is larger than the surge change. Notably, in R1 and R2, both surge and MSL change rates are positive, implying that their compound effect significantly contributes to increasing ESL. Several previous works have shown that MSL is rising along the Japanese coast (Senjyu et al 1999,  (1980-1999) and P2 (2000-2019) in Japan; (d) changes in the annual occurrence frequency of storm surge return period in P2 (2000-2019) relative to P1 (1980-1999). Sakurai 2006, Sasaki et al 2017), nevertheless, we report the compound effect of storm surge and MSL change for the first time.
When the annual mean latitude of peak surge location is analyzed over 1980-2019, there is a statistically significant northeastward change of 0.0549 [0.0607, 0.0533] degree per decade observed. Here, linear trend of µ from gamma distribution used to estimate the change rate as annual mean latitude of peak surge location did not follow normal distribution (Shapiro-Wilk test statistic: 0.9; p < 0.05; see figure S3 (b) in SI appendix). The northeastward shift of peak surge location likely leads to an increased surge levels in Tokyo (in R1), Yokohama (in R1), Nagoya (in R1), and Osaka (in R2) which are all located above the 34.5 • latitude belt. For example, the recorded number of surge annual maxima between the 34.5 • -35.5 • latitude belt is 1.7 times higher during P2 than during P1.
Following the assessment of changes in surge annual maxima and peak surge location, we also measure changes in the occurrence of different return period surge events ( figure 1(d)). It is particularly useful as it helps communicate whether environmental extremes are occurring more or less frequently relative to any chosen point in the past. First, we estimated the magnitude and annual occurrence frequency of a 2, 5, 10, 25 and 40 year storm surge event for 1980-1999 (see figure S4 in SI appendix). Second, we compare the annual occurrence frequency of similar return level events in P2 (2000-2019) with P1 (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999). Variations in different return periods' surge events mirror the changes in surge magnitude. A gradually increased frequency of small to mediumsized surge events has been observed in a recent 20 year period. For example, the ratio of the 25 year return period surge occurrence frequency in P2 to P1 is five. In Japan, the current storm surge warning criteria is based on the designed tidal level which is equivalent to a return period event (Hasegawa et al 2017). Therefore, our findings highlight the importance of regularly updating storm surge warning criteria under changing the return period of observed surge events.

Changes in TC meteorological variables influence storm surge
Here we discuss changes in TC meteorological conditions including track, intensity, size, and translation speed, which can significantly modulate storm surge climatology (Irish et al 2008, Rego and Li 2010, Sebastian et al 2014, Islam and Takagi 2020b. Besides the northeastward shift of peak surge location, there is evidence that TCs occurrence probability along the Pacific coast of the eastern (R1) and western (R2) Japan has increased substantially during P2 than in P1 (figure 2(a)). One remarkable feature is that the occurrence frequency of landfalling TC in R1 (nationwide) during P1 and P2 are 0.65 (1.85) and 1.15 (2.35) per year, respectively. It corroborates the findings of Kubota et al (2021), who detected a significant northeastward shift of TC landfall locations in Japan from 1977-2019. In contrast, TC occurrence frequency has remained almost stationary in R3.
Next, we analyzed the differences in TC landfall intensity (V max ) along with their six-hourly best track positions between P1 and P2. Overall, the frequency of strong TCs during the landfall time frame has increased in the last 40 years. For example, landfall TC activity with V max > 63 kt (∼equivalent to category 1 or higher category hurricane in Saffir-Simpson hurricane wind scale) has increased substantially in R1 and R2 from P1 to P2 ( figure 2(b)). In other words, the occurrence frequency of landfall TC intensity with V max > 63 kt has increased by 1.55 times. Similar increasing tendencies are also evident for other wind intensity categories (see figure S5 in SI appendix). These findings are consistent with Yamaguchi and Maeda (2020), who showed the increased frequency of strong TCs with a central pressure of less than 980 hPa which approached major Japanese cities (e.g. Tokyo, Nagoya, Osaka). It is also noticeable that TCs which made landfall in R3 have not changed their landfall intensity ( figure 2(b)). In addition to best track intensities, we examined recorded local maximum wind speed for each TC selected in this study (see figure S6 in SI appendix). It also depicts that TCs in P2 had caused stronger local winds (e.g. >50 kt) in R1, particularly than during P1. Figure 2(c) shows that the R2 region has been frequently exposed to large TCs (270 NM ⩽ R 30 ⩽ 430 NM), followed by R1 during landfall time in a recent 20 year period. Here, we adopted JMA's definition to describe the size of a TC (Large TC = 270 NM ⩽ R 30 ⩽ 430 NM; very large = R 30 > 430 NM; (JMA 2021a)) and provided the density map utilizing six-hourly best track positions along with TC size information. Overall, the frequency of large TC has increased more than fourfold from 1980 to 2019. TC-induced surge dependence over TC size is primarily due to sideway radiation and the dipole nature of the wind stress field. Both of them become less inhibitive for larger systems and, thus, result in greater surge levels (Nielsen 2009). However, it needs to be noted that a smaller but intense TC could also elevate surge levels as the pressure gradient tends to be steep, resulting in stronger winds near the eye of a TC (Islam and Takagi 2020a).
Along with analyzing TC track, intensity, and size, we also investigate the variability of TC translation speed. No apparent changes in frequency (no. of fastor slow-moving landfall TC) were observed during 1980-2019 (see figure S7 in SI appendix).
Although figures 2(a)-(c) demonstrate the statistical evidence of the changes in TC meteorological components individually that are likely responsible for modulating surge hazards in Japan, the amplitude of surge greatly depends on the combined effect of track, intensity, and size. Figure 2(d) illustrates the changes in the joint occurrence of strong (V max > 63 kt) and large (270 NM ⩽ R 30 ⩽ 430 NM) TCs over the past 40 years in Japan. It is noteworthy that Japan landfalling TCs are likely getting stronger and larger. This finding is consistent with Lin and Chou (2018), who showed that TCs traversing the Philippines are getting stronger and larger during pre-landfall time (0-48 h) based on a 37 year  JMA best track dataset. Observational changes in surge annual maxima in R1 and R2, in general, reflect the combined effect of shifting TC landfall location, increasing intensity and size of TC. More explicitly, the integration of storm size and wind strength over the footprint of TCs provides a bulk amount of energy/momentum transferred from the storm to the water column and, thus, the functional dependence of the surge level on the velocity and storm radius. The environmental conditions include the sea surface temperature, vertical wind shear, and relative humidity at 500 hPa, atmospheric circulation conditions, and anthropogenic climate forcings may be responsible for creating favorable conditions for such kind of strong and large TC development in P2 than in P1 (Liu andChan 2020, Yamaguchi andMaeda 2020); however, a detailed analysis of these processes is beyond the scope of this study.
The possible association between the change in surge annual maxima and TC meteorological condition has been discussed above from an overall perspective. Quantifying individual and joint contributions of TC meteorological variables is also important to understand which variables are likely influential to the observed surge changes. Therefore, we applied the SSHPI to demonstrate it quantitatively. Note that the SSHPI is correlated to the observed  (1980-1999; n = 37). Six-hourly JMA best track positions along with TC wind intensity (Vmax) and size (R30) information are utilized to provide density maps. peak surge height and explains ∼57% of the observed variance ( figure 3(a)). These statistics are comparable to storm surge estimation obtained from full physical numerical surge models (e.g. Dawson 2014, Marsooli and. Thus, SSHPI is a reasonably good proxy for TC-induced peak surge threat and could be utilized to quantify the contributions of TC meteorological components to hazardous surge height. To estimate the possible contributions of observed TC meteorological variables over the changes in surge annual maxima, we separately changed two SSHPI components (V max or R 30 ) or a set of components (V max and R 30 ) to their P2 (2000-2019) averaged values and kept the other components at their P1 (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999) averages. Finally, we estimated each change that affected total surge change by comparing P2 and P1 index values (see table S1 in SI appendix). In R1, the changes in TC intensity and size contributed to the changes in surge annual maxima by 41% and 21%, respectively and their joint contribution may be responsible for as much as 71% ( figure 3(b)). This joint influence statistic gets smaller in R2 (23%) and it has no significant role in R3. It should be noted that compared to R1, the influence of TC intensity and size over the changes in surge level in R2 is limited. This may be because other influencing variables, such as wave set-up and inverse barometer effect are responsible which were not considered in the SSHPI. Nevertheless, figure 3(b), in general, mimics TC meteorological conditions illustrated in figure 2. Thus, it is reasonable to believe that ESL differences between the two 20 year periods studied here cannot be explained by MSL rise alone, rather the northeastward shifting of TC landfall location, as well as intensifying and widening of TC might have altered the likelihood of ESL including surge extremes. These findings are not consistent with previous global/regional ESL trend studies (e.g. Woodworth and Blackman 2004, Church et al 2006, Menéndez and Woodworth 2010 in which surge trends are masked by MSL change. While MSL change is independent of TC activity, storm surge hazard is sensitive to local scale change in TC meteorological variables. Therefore, it is not implausible to observe significant variability of surge hazard on a local scale. In addition to the long-term changes in MSL, we also emphasize that coastal adaptation and mitigation strategies for any future ESL changes should consider middle-short term variability of TC frequencies, characteristics, and resultant surge hazards.  (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999) while maintaining the other variables value at their averages from P1.

Conclusion
Here we have addressed notable positive changes in surge annual maxima in eastern and western Japan during 1980-2019. It is likely due to the northeastward shifting of TC landfall location, as well as intensifying and widening of TC over the recent 20 years (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019). Contrary to what previous studies suggest, changes in surge levels were comparable to those observed for MSL rise, leading to the overall change in ESLs over the last 40 years. The Japanese coast will likely experience the increase of extreme storm surge events if the current trend of TC intensification and widening continues.
This study statistically quantifies the contribution of TC meteorological variables including wind speed, size, and track of TCs in changing surge hazards; nonetheless, other factors such as air pressure and wave set-up can also modulate surge intensity. The physical meaning of these influential factors remains unknown. Therefore, the surge trend attribution will likely require further research including model evaluation on the contributions of TC meteorological variables, anthropogenic forcings, and interdecadal variations using numerical simulations. Furthermore, expanding the scope of this current study to include storm surges caused by non-landfalling TCs that passed near the Japanese coastlines and extra-TCs, could serve as a potential direction for future research.

Data availability statement
The data cannot be made publicly available upon publication because they are owned by a third party and the terms of use prevent public distribution. The data that support the findings of this study are available upon reasonable request from the authors.