Hybrid methods combining atmospheric reanalysis data and a parametric typhoon model to hindcast storm surges in Tokyo Bay

Reanalysis data and a parametric typhoon model formula are often used to prepare wind and pressure fields for storm surge hindcasting. However, their optimum selection and usage have not been well established. To enhance the accuracy of wind and pressure fields, two hybrid methods were proposed by applying a parametric typhoon model of Mitsuta–Fujii’s formula, which is determined from the typhoon center to a certain radius of Rb, and then switching to the reanalysis data of ERA-Interim in the outer region through the interpolated transition bandwidth of Wb. In hybrid model I, Rb and Wb were fixed because the time-varying radius of the maximum wind speed was determined by the typhoon formula. In hybrid model II, these parameters were determined to minimize the mean difference between the reanalysis data and the fields obtained by typhoon formula in the transition band at each time step of the typhoon track. The hindcasting of eight significant typhoon events approaching Tokyo Bay was performed. This validated the proposed methods in comparison with the observed storm surge anomalies. Both models performed satisfactorily. Hybrid model II was found to be superior in terms of the balance of accuracy and preparation cost.


Results
By introducing the radius of R b and the transition bandwidth of W b , a generalized hybrid method is proposed where the M-F model is applied in the 0 ≤ r ≤ R b region and switching to ERA-I in the r ≥ R b + W b region (outer region of the transition band) while interpolating the two models in the transition band of R b < r < R b + W b . In Section 2.2, this generalized idea (hybrid model II) is introduced with a method to determine the R b and W b parameter values. Before that, as a special case, hybrid model I is introduced in Section 2.1 by applying the M-F model in the region between the typhoon center (r = 0) and the applicable boundary, r = r max , and switching to www.nature.com/scientificreports www.nature.com/scientificreports/ ERA-I in the region where r ≥ 2r max . ERA-I and M-F model were smoothly interpolated in the transition band of r max < r < 2r max .
Hybrid model I. The direct method of determining the wind and atmospheric pressure field is given by: where F B is the resultant wind velocity component in x or y direction, or the atmospheric pressure, F M−F and F ERA are the wind velocity or atmospheric pressure from M-F model and ERA-I, respectively, r is the distance between the target location and typhoon center, and r max is given by equation (6). Compared with the method of Shao et al. 27 , the proposed hybrid model can be easily used, and the weighting coefficient is explicitly explained. Comparing the wind velocity distributions between the ERA-I and hybrid model I (see Supplementary  Fig. S1), ERA-I and the M-F model dominated outside and inside the typhoon region, respectively, and the transition between these two models was continuous.
Comparisons between the computed (ERA-I, M-F model, and hybrid model I) and observed (JODC) wind speeds are shown in Supplementary Fig. S2. Their coefficients of determination are presented in Table 1. It can be seen that wind speeds computed by hybrid model I were more consistent with those of the observed (JODC) values during Typhoon 8506 and Typhoon 1115 than those computed by the M-F model or ERA-I.
Hybrid model II. The second hybrid model is based on the analysis of the difference between ERA-I and the M-F model. As shown in Section 2.1, the accuracy of the M-F model is higher than that of ERA-I around the typhoon center while the accuracy of ERA-I is higher than that of the M-F model away from the typhoon center. Thus, an optimum switching method between the M-F model in the central part and ERA-I in the outer typhoon region was obtained by interpolating the two data in the transition band at each time step for the typhoon track. At each time step, both the distance r = R b from the typhoon center to the inner side of the transition band and the bandwidth W b were investigated from r = 0 until r ≤ L max (L max is the maximum searching length) to minimize the mean value of the differences between ERA-I and M-F model across the computational grids within the transition band. For simplicity, a variable of the searching increment, Δr, was introduced where R b = (m − 1)Δr (m is an integer number and 1 ≤ m ≤ n) and the maximum integer number of n was determined to satisfy (n − 1) Fig. 2). The formula for hybrid model II is given by: The increment value of Δr was set as 5 km and the maximum searching length of L max was set as 1000 km, considering the balance between the accuracy and computational cost. Supplementary Fig. S3 shows the time series of the optimum R b to r max ratio for ten cases of W b to obtain the optimum bandwidth. The ratios were always 0 when the distances from the typhoon center to the computational nodes were greater than 1000 km, e.g., in the period until 18:00 on June 28 th . For all ten cases, the ratios were generally less than 2, which is consistent with the applicable radius, r, range for the M-F model (0 ≤ r ≅ r max ) as suggested by hybrid model I. To further evaluate the performance of hybrid models I and II, eight selected typhoon cases were analyzed, and their corresponding root mean square errors (RMSE) were calculated by: where F hb−I,i is the wind velocity value or atmospheric pressure at the i th time step during a certain typhoon computed by hybrid model I, F JODC,i , is the corresponding measured wind speed or atmospheric pressure, and N is the total number of typhoon time steps. Supplementary Table S1 presents comparisons of the RMSE values among ERA-I and hybrid models I and II with different bandwidths. Compared with the original ERA-I RMSE values, both hybrid models were found to improve accuracy. However, the results from other cases showed that hybrid model II performed better with varying optimal bandwidths, compared with hybrid model I. Hybrid model I was considered a special case of hybrid model II where the bandwidth and band distance from the typhoon center were determined based on ERA-I and M-F model. Although the computational cost for hybrid model II was slightly higher, its performance was better and thus, would be a better choice for enhancing the accuracy of storm surge computation. The hybrid model II concept was then applied to the atmospheric pressure fields, which is also vital in enhancing the storm surge computation accuracy. Comparisons of the time series of atmospheric pressures at Station Tokyo and Station Yokohama for Typhoon 8506 and Typhoon 1115 between ERA-I, Myers formula, hybrid model II, and observed data are shown in Supplementary Fig. S5. The blended atmospheric pressures obtained by hybrid model II were more consistent with those of the measured values. Thus, hybrid model II is confirmed to be superior.

Application to hindcast storm surge.
To validate the performance of the proposed hybrid models for wind and atmospheric pressure fields, a hindcasting of storm surges in Tokyo Bay was performed using the finite volume community ocean model (FVCOM) forced by the wind and pressure data of ERA-I, M-F model with Myers formula, and the two hybrid models. To further consider the effect of possible water volume exchange Searching band areas in hybrid model II (the first parameter is the increase in the searching distance Δr, the second parameter is the maximum searching radius L max , the third parameter is the bandwidth W b , which is the equal for all searching bands, and the fourth parameter is the distance from the typhoon center to the inner side of the transitional annulus band area L i (i = 1, 2, …, n). R b equals L m if the m th annulus is the band area within which the difference of ERA-I data and the M-F model is minimum. The searching begins from the first circle (black color) i = 1 with radius W b , then the blue color annulus i = 2 with internal radius L 2 = Δr, etc. The investigation stops when L n + W b = 1000 km. It is clear that in the proposed searching method, part of the neighboring band area is overlaps because the searching band area gradually expands with Δr). (2019) 9:12222 | https://doi.org/10.1038/s41598-019-48728-7 www.nature.com/scientificreports www.nature.com/scientificreports/ between the ocean from a far place and the inner side of the bay, and also to reduce the boundary influence, a wide area was selected as the computational domain (20°N-60°N, 120°E-160°E) as shown in Fig. 3. A total of 34,255 computational nodes and 61,905 triangular elements were used.
Supplementary Table S2 summarizes the information on the storm surge observation stations 38,39 , and Fig. 1a shows the station locations in Tokyo Bay. To compare the performance of the hybrid wind data, the same set of parameters and boundary conditions were adopted except for the wind and atmospheric pressure data. The computation started with the still water condition with a mean sea level and zero velocities across the domain. For the wind data, ERA-I, M-F model, hybrid model I, and hybrid model II were used. For the atmospheric pressure data, Myers formula, ERA-I, and hybrid model II were applied. Supplementary Fig. S6 shows the comparisons of the storm surge anomalies between the hindcasting results and observed data at six stations in Tokyo Bay for Typhoon 8506. Before the storm surge peaked at 5:00 on July 1 st , 1985, the observed sea level anomalies at six stations were generally higher than the mean sea level (=0). ERA-I reproduced the observed tendency better than the M-F model with Myers formula. This indicates the significant influence of the atmospheric pressure fields on the reproducibility of sea level anomalies before the typhoon approach. However, ERA-I data substantially underpredicted the storm surge peak at all six stations, which implies the insufficient spatial resolution of ERA-I data to reproduce the surge peak. Comparing all the model results, hybrid model II (for both wind and atmospheric pressure fields) performed better in reproducing the storm surge anomalies during the typhoon approach, the surge peak, and the following oscillations induced by bay-scale resonance. Supplementary Fig. S7 shows comparisons of the storm surge anomalies between the computed values using six models and the observed values for Typhoon 1115 at four stations in Tokyo Bay. M-F model could not reproduce the storm surge anomalies well before the peak, and the peak anomalies in Station Tokyo Harumi, Station Chiba, and Station Yokosuka were significantly lower than the observations. Similarly, ERA-I could not reproduce the peak anomalies at four stations. Among the six storm surge models, only hybrid model II (for both wind and atmospheric pressure fields) exhibited better results in reproducing storm surge anomalies.
To further evaluate the accuracy among the six models, the RMSE (defined by equation (3)) between the computed and observed values are presented in Table 2 for Typhoon 8506 and Typhoon 1115. The following are the periods of comparisons: from 6:00 on June 29 th to 9:00 on July 2 nd , 1985, for Typhoon 8506, and 3:00 on September 20 th to 23:00 on September 23 rd , 2011, for Typhoon 1115. According to Table 2, it can be concluded that the wind and atmospheric pressure fields obtained by hybrid model II should be used to hindcast storm surges.

Discussion and conclusions
Parametric typhoon models and reanalysis datasets are widely used for storm surge simulation. However, only a few studies seek their optimum combination to enhance the accuracy of storm surge computation. In this study, two hybrid models were proposed for wind fields, where the M-F model was applied between the typhoon center (r = 0) and a certain radius (r = R b ), and switched to ERA-I when the radius is greater than R b + W b through the linearly interpolated transition band with the width of W b . For hybrid model I, both R b and W b were fixed as the time-varying radius of r max was determined by the M-F model at a distance from the typhoon center where the maximum wind speed occurred. Conversely, in hybrid model II, the optimum combination of R b and W b was determined to minimize the mean difference between the M-F model and ERA-I within the transition band. www.nature.com/scientificreports www.nature.com/scientificreports/ Thus, hybrid model II was considered to be the generalization of hybrid model I. The wind fields of the eight historical typhoon cases approaching Tokyo Bay were compared between ERA-I, the M-F model, hybrid model I, and hybrid model II. The accuracy of the typhoon wind fields calculated by ERA-I, the M-F model, hybrid model I, and hybrid model II were verified by comparison with the observed wind field provided by JODC. Results showed that while both hybrid models performed better than ERA-I and the M-F model, the accuracy of hybrid model II was higher than that of hybrid model I. The atmospheric pressure fields were also modified using hybrid model II.
Using FVCOM, hindcasting of storm surges in Tokyo Bay was performed by forcing six combinations of modeled wind and atmospheric pressure fields. The results showed that the modification of both atmospheric pressure and wind fields significantly improved the accuracy of the storm surge anomalies. Both proposed hybrid models performed better than the computations using only ERA-I or the M-F model. Hybrid model II can be easily tuned for each storm surge cases using only ERA-I and the M-F model data, which significantly improved the accuracy of storm surge hindcasting.
Forecasting the storm surge a short time period ahead of the typhoon's arrival was out of the scope of this paper. However, the proposed method may be applied to the storm surge prediction under future climate conditions. The future atmospheric fields could be prepared using a general circulation model (GCM) with future climate scenarios, e.g., scenarios discussed in Intergovernmental Panel on Climate Change (IPCC) reports. Then, these atmospheric fields can be used in place of ERA-I reanalysis data for hybrid model I and II combining with parametric typhoon model.

Typhoon cases and observation stations. The locations of the meteorological observation stations in
Tokyo Bay and its surrounding areas are shown in Fig. 1a. The stations recorded the tide levels, including the storm surge anomaly, and wind speed in a one-hour time interval. The detailed information, including the station names, longitudes, and latitudes, is summarized in Supplementary Tables S2 and S3. Among historical typhoons and storm surges that occurred in Tokyo Bay between 1951 and 2017, 15 typhoon cases with large storm surge anomalies (greater than 0.8 m) were screened out and eight typhoons (see Fig. 1b) were selected as summarized in Supplementary Table S4.
Wind data analysis. Several reanalysis datasets are for wind and pressure fields, including the ERA-I, NCEP-DOE reanalysis II, and parametric typhoon models. A preliminary investigation was performed on the data consistency by a comparison with observed data from eight typhoon cases obtained in Section 4.1. ERA-I and the M-F model were used because both of them were more consistent with the measured data for the eight typhoon cases in Tokyo Bay during the preliminary study. Moreover, this method is applicable to any combination of reanalysis datasets and typhoon formulas.
ERA-Interim reanalysis data. ERA-Interim (ERA-I) is a global atmospheric reanalysis production provided by the European Center for Medium-Range Weather Forecasts (ECMWF), which started in 1979 and has since been continuously updated. In this study, ERA-I data of 6-hour interval wind speeds at a 10-m height from the mean sea level and the sea surface level atmospheric pressure were used. One-hour interval time series datasets for ERA-I were also created by interpolation. To obtain the most accurate ERA-I at the target stations, a 0.125° grid in the longitude and latitude was used. Figure 4 shows a comparison of ERA-I and the measured wind data provided by Japan Oceanographic Data Center 40    With the typhoon center far from the stations, the wind speeds were fairly consistent with the observed data with the corresponding approaching and departing distances between the typhoon center and Station Tokyo, Station Honmoku, and Station Shionomisaki ranging from 3.94 r max to 1.84 r max (1985/7/1/0:00-1985/7/1/10:00), 3.09 r max to 2.23 r max (1985/6/30/18:00-1985/7/1/10:00), and 2.15 r max to 2.44 r max (1985/6/29/12:00-1985/7/1/12:00), respectively. With the aforementioned distances, large discrepancies between the observed data and ERA-I were noted.
After the ERA-I analysis for the eight typhoon events and the JODC observation data at eight stations, the applicable ranges for ERA-I are summarized in Table 3. The critical value, D 1 , is defined as the applicable distance limit from the typhoon center. It was proposed to identify the applicable range of ERA-I at a target location during the typhoon approach. ERA-I is applicable at a location when the distance from the typhoon center to the location is greater than D 1 . A similar explanation can be applied to D 2 , which is defined for the typhoon departure condition. The distances D 1 and D 2 consist of the applicable boundaries for ERA-I.
The first step in determining D 1 and D 2 values is the plotting of the time-series ERA-I data and JODC observed data at each station during the selected typhoon events. Instances with large discrepancies are then noted. The distance from the typhoon center to the station at the noted instance is finally calculated as the D 1 value for this typhoon. The distance of D 2 is similarly determined under the departing condition. As shown in Table 3, different stations experienced the whole processes of the eight selected typhoon events. Statistics at the eight stations under different typhoon events resulted in different D 1 values. It can be seen that D 1 and D 2 values are approximately 2 r max . Hence, to simplify the following analysis, D 1 and D 2 are set as 2 r max .
Parametric typhoon model description. The Myers formula based on the exponential distribution of the atmospheric pressure field is given by: where P(r) is the pressure at a radial distance r from the typhoon center, P c (hPa) is the typhoon central pressure, P 0 (=1013.25 hPa) is the ambient or environmental pressure, r is the distance from the computational mesh node to the typhoon center, and r max (km) is the maximum wind speed radius. After reviewing similar studies 17,41 , the M-F model was selected to compute the wind field as presented in equation (5) and the estimated pressure by equation (4) was applied to the M-F model. www.nature.com/scientificreports www.nature.com/scientificreports/ w1 max T U w is the total wind vector, U w1 is the moving component, U w2 is the wind vector induced by the rotating component, P(r) is the pressure field calculated by Myers formula, C 1 and C 2 are dimensionless coefficients ranging from 0.6 to 0.75, f is the Coriolis parameter, r is the distance from the typhoon center, ρ a is the atmospheric density, and V T is the typhoon forward speed obtained from the best track data 29 provided by the JMA including the typhoon center location and central pressure. The time-varying radius of the maximum wind speed, r max , was determined as a function of central pressure, P c , following the empirical formula 17  Applying this model to the eight typhoon cases, the estimated wind speeds were compared with the observed data for Typhoon 8506 and Typhoon 1115 cases as shown in Supplementary Fig. S2. In general, for large wind speeds, the estimated values were consistent with that of the measured data. Following Section 4.2.1, the applicable boundary distances of D 3 and D 4 were determined for the M-F model from the typhoon center to the target location when the typhoon is respectively approaching and departing, as shown in Table 3 (D 3 and D 4 are approximately equal to r max ). Thus, M-F model is applicable in the area between the typhoon center (r = 0) and the boundary (r = r max ).  Table 3. Critical values of D 1 -D 4 at eight stations under eight selected Typhoons.