Bias correction of simulated storm surge height considering coastline complexity

In this study, we propose a new approach for model validation that can be applied to the projection of possible future storm surge heights (SSHs) on the regional scale. First, this study conducts a series of SSHs for the southeastern coast of the Korean Peninsula (KP) by six typhoons that produced SSHs over 1.0 m since 1979 and identifies the bias between simulated and observed SSHs. Next, formulas for the bias correction using a geographic parameter, in particular the coastline complexity factors, are drawn and validated. Finally, the effect of the proposed bias correction on projection of future SSHs is examined by performing simple tests to consider only central pressure drops to reflect the impact of climate change. It can be seen that the bias correction method considering the coastline complexity can improve the model’s accuracy by 14% to 23% and prevent potential overestimation by up to 20% of the maximum SSHs considering climate change effect on the southeastern coast of the KP.


INTRODUCTION
According to previous studies, storm surge heights (SSHs) are expected to become more severe in the future climate condition because climate change is projected to increase sea surface temperature and consequently results in intensified tropical cyclones (e.g. Emanuel, 2005;Knutson et al., 2010;Mori and Takemi, 2016). The Korean Peninsula (KP) in North-East Asia is no exception (e.g. Kim et al., 2014;Oh and Moon, 2013). Therefore, the projection of potential maximum SSHs and the identification of areas vulnerable to storm surge on regional or national scales are important for policy making to prioritize where to take action and to find appropriate site-specific solutions against the SSHs in the future climate condition along the coastlines (Vousdoukas et al., 2016;Lapidez et al., 2015).
In order to project the potential maximum SSHs on a regional or national scale in the KP, it is necessary to perform a series of simulations of SSHs by various typhoons Correspondence to: Jung-A Yang, Graduate School of Engineering, Kyoto University, Nishikyo-ku, Kyoto 615-8530, Japan. E-mail: yangja.1985 @gmail.com that have different properties such as tracks, magnitudes and moving speeds including consideration of climate change effects. Prior to evaluating these SSHs using a numerical surge model, it is necessary to validate the model. Many researchers have validated numerical models by adjusting the values of the physical parameters in the model to obtain simulated SSH values close to the observed values for a particular typhoon at a particular location/site (e.g. Shibutani et al., 2015;Nakajo et al., 2015). However, it is not generalized for other locations, events and different model configurations. This is because hindcasting a specific typhoon-induced storm surge has uncertainties depending on topographic features, typhoon characteristics and model configurations. Moreover, the long-term impact assessment of storm surge requires a large number of ensembles. Therefore the spatial resolution of computational mesh needs to be coarse to reduce computational cost. It is necessary to implement general bias correction of numerical results considering both spatial resolution of mesh size and coastal features for impact assessment of climate change on storm surge study.
Therefore, this study proposes a new approach for model validation that can be applied to the projection of future SSHs on the regional scale. We first simulated SSHs on a regional scale in the KP for various historical typhoons and examined the performance of the model. Next, the bias between the simulated and observed SSHs was modified to take into account a geographic parameter, in particular the coastline complexity of each individual location/site. Finally, the effect of the bias correction on projections of the potential maximum SSH was examined through a simple test.

Selection of target typhoons
According to the National Typhoon Center (2016), the KP was impacted directly and indirectly by 345 typhoons from 1904 to 2015, on average three typhoons per year. As it is too expensive to analyze all storm surge events resulting from all typhoons that affected the KP historically, this study considered six typhoons that produced SSHs over 1.0 m in the KP.
Forty-seven tidal stations are located along the coast of the KP (shown by red circles in Figure 1) and have been operated by the Korea Hydrographic and Oceanographic Agency since 1952. The six typhoons were selected based on the highest SSHs observed at all tidal stations in the KP from 1979 to 2013. The observed SSH is in general obtained by subtracting astronomical tide from observed tidal data. We extracted the astronomical tide from the observed tidal data using harmonic analysis. After calculating the observed SSHs, it turned out that all observed SSHs higher than 1.0 m occurred in the southeastern coast of the KP. These six typhoons were named Typhoons Maemi (2003), Rusa (2002), Bolaven (2012), Nari (2007), Megi (2004) and Sanba (2012) (see Table I and Figure 1). We selected these six major typhoons as the target typhoons

Numerical model and model setup
The simulated SSHs by the target typhoons were calculated using a coupled surge, wave and tide model (SuWAT) developed by Kim et al. (2008Kim et al. ( , 2015. The target area has complex topography but gentle bathymetry gradients. Therefore, we considered only the surge module to reduce the computational cost. For atmospheric forcing, an empirical typhoon model embedded in the SuWAT model was used. To generate wind and pressure fields, the typhoon parameters were read from the National Institute of Informatics (2014), except the radius of maximum wind (R max ) which is given by an observed empirical formula (Yasuda et al., 2010). R max is an indicator of the size of a typhoon and has different values for each typhoon.
We carried out a series of storm surge simulations, one for each target typhoon. Each simulation was run with 48 hours spin-up prior to the calculation of the typhoon event to bring stability in the sea surface level over the entire computational domain. Four-domain nesting was used with a 1:3 ratio of spatial grid size from coarse grid to fine grid (shown by the dashed rectangles in Figure 1). The outermost domain (D1) is set in the region of 20°-42°N, 117°-132°E that almost covers the formation and dissipation of target typhoons and the innermost domain (D4) is the region 33.5°-35.75°N, 127°-130°E which we define as the southeastern area of the KP and has a spatial resolution of 0.008333° (about 760 m). The bathymetry data in D1 were given from the General Bathymetric Chart of the Oceans (GEBCO, 2003). For the remaining nested domains, the bathymetry data were derived from the Digital 30 s Gridded Bathymetric Data of the Korea Marginal Seas dataset (Seo, 2008). The time step of storm surge computation was 300 s. The inundation was not considered in this study.

Bias correction method considering coastline complexity
The simulated SSHs over this region have some disagreement with observed SSHs of the target typhoon, which we attributed to simulation model bias. In this section, we present a formula for bias correction using a geographic parameter, in particular the coastline complexity factors. First, we drew a relational expression between bias and coastline complexity factors. The observed SSH (ζ obs ) may be expressed by the sum of the calculated SSH (ζ cal ) To express Δζ as an equation, we look to an empirical model for SSHs. The Japan Meteorological Agency used the following empirical formula for predicting storm surge height before a numerical model was developed (Unoki and Isozaki, 1966): where ζ predict [m] is the predicted SSH from the empirical model, Δp [hPa] is the difference between environmental atmospheric pressure (1010 hPa) and lowest atmospheric pressure at the typhoon location, U [m/s] is the maximum wind speed, θ [°] is the angle between predominant wind direction that would cause the highest SSH at a location and the wind direction at the time of wind speed U, and a [m/hPa] and b [s 2 /m] are the tuning coefficients determined for the location based on historical observation data. Generally, regardless of location, 0.0099 is used as the value of the coefficient a, but the coefficient b strongly depends on regional features such as characteristics of the coastline, the shape of a bay and local bathymetry (see e.g. Miyazaki et al., 1961). That is, the coefficient b is an inherent factor reflecting the characteristic of terrain and can be expressed where ζ pressure [m] and ζ wind [m] are sea level rise by atmospheric pressure drop and by wind, respectively. For regional correction of the empirical model, the geographical parameter b should be replaced with (b + Δb).
We assume ζ predict is comparable to ζ cal : If we can apply the above method to the numerical result, we can express the calculation error as the equation: Here, we call Δb a terrain correction factor and assume Δb is the difference between b obs and b cal , calculated by substituting the ζ obs and ζ cal to ζ predict in Equation (3), respectively: However, Equation (7) cannot be used to project the future SSH, because it needs the observed SSH to estimate Δb. So, we parameterized Δb as a function of the coastline complexity factor F D , where, f(F D ) [s 2 /m] is a tuning function of the F D which is a dimensionless parameter. The cos θ term in Equations (6) and (8) is taken to be one because we are concerned only with the maximum wind-induced surge in this paper (i.e. where wind direction is perpendicular to bay direction). Next, the F D at each location in the target area was calculated using the concept of fractal dimension and then Δb is expressed as a function of F D by deriving regression equations between them. The fractal dimension was introduced by Mandelbrot (1967) and has been used as a measure of the complexity of objects. Many researchers have applied the fractal dimension to assess complexities of coastlines in different countries: Britain (e.g. Philip, 1994), USA (e.g. Carr and Benzer, 1991) and China (e.g. Zhu et al., 2004). In this study, the box-counting method was adopted to determine the coastline complexities in the target area. As base data of coastline shape, images of coastline maps based on version 2.3.2 of the Global Self-consistent, Hierarchical, High-resolution Shoreline Database (Wessel and Smith, 1996) were used for analysis.
Finally, the verification of this bias correction method and the estimation of the effect on the projection of potential maximum SSHs were conducted.

Hindcast of storm surge
In Figure 2, the simulated maximum SSHs for historical events are compared with the observed maximum SSHs available from the nine tidal stations shown in Figure 3. When the distance between the center of the typhoon and the tidal station is small and the tidal station is located to the right side of the typhoon, the SSHs tend to be overestimated (refer to Figure S1). This is especially pronounced among tidal stations that have complex coastlines. However, when the distance between them is large, the SSHs tend to be underestimated regardless of relative position between them. The more complex the coastline, the larger the error range. So, it is necessary to consider topographical features to correct the bias. Furthermore, the overestimation (underestimation) appeared to be related to values of SSH higher (lower) than 1.0 m. The model performance at each tidal station was evaluated using the root mean square error (RMSE) with values from 0.13 m at Masan station to 0.64 m at Tongyeong. The RMSE values tend to be larger for tidal stations located along more complex coastlines. The RMSE for overall model performance is 0.43 m. The reason why the RMSE is large may be because we calculated the storm surge heights for each target typhoon using the same values for all physical parameters in the numerical model.
So far, many researchers have conducted validations of numerical models by hindcasting historical SSHs using historical typhoon data and observation surge height data. However, most researchers verified the numerical models not by correcting calculation errors but by adjusting values of physical parameters to obtain SSH close to the observed one. For example, Shibutani et al. (2015) verified the numerical model by hindcasting SSH at Ise Bay, Japan, during the typhoon Vera (5915) event. In their study, they adjusted the value of reduction factors for wind speed at 10 m relative to the typhoon's movement and gradient wind speed. This parameter is used to reduce the influence of the wind velocity on the right side of the typhoon path when calculating the wind field by the parametrical wind model. Therefore, as the path of the typhoon changes even for the same region, the proper value of that parameter can vary. In addition, when estimating SSHs at multiple points by a sin- gle typhoon, the optimum value of that parameter can be different at each point. In another case, Nakajo et al. (2015) modulated the value of R max to confirm the reproducibility of the numerical model by hindcasting SSH that occurred in Tokyo Bay, Japan, during typhoon Roke (1115). Thus, previous approaches to model validation are suited for the estimation of SSH at a particular point for a particular typhoon, but these approaches are not suitable for the estimation of SSH at multiple points for various typhoons.
Henceforth, we examined a new bias correction method considering coastline complexities that can be applied to the study on the projection of possible future change in SSHs on the regional scale. The results for four of the six typhoons were used to derive the bias correction coefficients and the results of the remaining two were used for validation.

Tuning of the bias
For the bias correction, first, the fractal dimension F D was estimated on each grid in D4 that included the nine tidal stations. The F D at the nine tidal stations ranged from 1.020 to 1.268 (see in Figure 3). Since larger values of F D signify more complex coastlines, among the nine tidal stations the Tongyoung station has the highest coastline complexity.
Next, we derived regression formulas between the F D of the nine tidal stations and the coefficient Δb obtained from storm surge simulations for Typhoons Meami, Megi, Nari, and Sanba. In Figure 2, the number of observation data above 1.0 m is very limited. In that range, the error range tends to also be large. These four typhoons were chosen to draw the regression equation because they have high SSHs and large errors, and the resulting bias correction performs well. A linear regression equation (Equation (9), abbreviated as LRE) and a cubic regression equation (Equation (10), abbreviated as CRE) were drawn as follows: By conducting the bias correction based on LRE and CRE, the RMSE of the original results apparently decreased from 0.40 m to 0.33 m and to 0.28 m, respectively (Figures 4ac). The bias correction by LRE and CRE improves the model's performance by 18% and 30%, respectively.
Finally, these regression formulas were validated using results of Rusa and Bolaven events. The RMSE of results without bias correction was reduced from 0.51 m to 0.45 m (11%) and 0.44 m (14%) by applying the LRE and CRE for bias correction, respectively (Figures 4d-f). However, the influence of the higher order regression equations on model performance is not large. By observation point, in all tidal stations except Yeosu the RMSE of results modified by LRE (CRE) decreased by -0.04 m to -0.16 m (by -0.01 m to -0.18 m) and the amount of decrease tended to be larger at tidal stations with simple coastlines. At Yeosu station, the RMSE value was increased (decreased) by 0.20 m (-0.18 m) when LRE (CRE) was applied for bias correction.
In general, the RMSE of results without bias correction was reduced from 0.43 m to 0.37 m (14%) and 0.33 m (23%) by applying the LRE and CRE for bias correction, respectively. The CRE is more effective than LRE. It can be seen that the bias correction method considering the coastline complexity can improve the model's accuracy.
To confirm the importance of considering terrain complexity for bias correction, we compared the results with bias correction for all available data using the proposed method and using a conventional method using the average difference between observed and simulated SSHs at various SSH levels. The RMSE of each result is 0.33 and 0.32, respectively (refer to Figure S2). The small surge magnitude mainly occurred due to pressure induced surge, therefore, the bias is rather linear. Both bias correction methods can give reasonable improvement for the small surge region. The RMSE for the small surge region is 0.28 and 0.23, respectively. However, the significant surges with height greater than 1 m mainly occurred due to wind induced surges. In these cases, the relational position between typhoon track and locations of interest, terrainspecific characteristics such as coastline complexity, the angle of wind to target points and bathymetric features result in different characteristics of peak surge as shown in Figures 1 and 4 in the manuscript. The RMSE for the large surge regions is 0.43 and 0.49, respectively. Additionally, we attempted a third bias correction method using the average ratio between the observed and the simulated SSHs. The RMSE of this method is 0.48 (refer to Figure S3c). This method is not effective for bias correction. Our proposed bias correction method considering coastline complexity gives a smaller RMSE value in the large surge region than the two alternatives. However, further refinements of the proposed bias correction method considering not only coastal complexity but also other characteristics such as depth and sea bottom slope remain the next challenge.

Evaluation of the effect of the bias correction on projection of future storm surge
To assess the impact of climate change on storm surge, several approaches exist such as dynamic approaches using the output from a general circulation model or regional climate model (e.g. Yasuda et al., 2014), statistical approaches based on a stochastic typhoon model (e.g. Nakajo et al., 2013) and worst case approaches based on extreme scenarios (e.g. Shibutani et al., 2015). But, there is no quantitative consensus on projections of typhoon intensity in the litera- Here we performed simple sensitivity tests to select an arbitrary number for central pressure drops to reflect the impact of climate change. We designed three future scenarios in which the central pressure drops for the target typhoons were assumed to be 10 hPa (future scenario, abbreviated as FS 1), 20 hPa (FS 2) and 30 hPa (FS 3). All other physical parameters are taken to be the same as in the present scenarios.
The simulated maximum SSHs for all target typhoons under these simple future scenarios showed similar trends. Due to limited space, only the potential maximum SSHs by latitude for the Typhoon Rusa event are shown here in Figure 5 (a). It can be seen that as the central pressure drop becomes larger, the potential maximum SSHs along the coast becomes higher as expected. The potential maximum SSHs modified by LRE and CRE are plotted in Figures  5 (b) and (c), respectively. By applying the bias correction by LRE (CRE), the largest potential maximum SSHs decreased about 22% (12%) from 2.59 m to 2.01 m (2.27 m) under the present condition, by 21% (13%) from 3.21 m to 2.52 m (2.80 m) under the FS 1, by 21% (13%) from 3.85 m to 3.04 m (3.34 m) under the FS 2 and by 20% (13%) from 4.52 m to 3.61 m (3.92 m) under the FS 3. From these results, we conclude that the bias correction considering the coastline complexity could prevent overestimation of the potential maximum SSHs on the southeastern coast of the KP.

CONCLUSIONS
This study proposed a new approach for model validation that can be applied to the projection of potential maximum SSHs on the regional scale.
First, a series of storm surge simulations were conducted for the southeastern coast of the KP by target typhoons. The simulated SSHs over this region had some disagreement with observed SSHs, which we attributed to simulation model bias. Next, formulas for the model bias correction using the coastline complexity factors were drawn and validated. The simulated results for four of the six typhoons were used to derive the bias correction formulas. The proposed bias correction formulas improved the model's performance by 18% to 30%. In validation, the simulated results of the remaining two were used and the RMSE of results without bias correction was reduced by 11% to 14% by applying the bias correction method. Overall, the RMSE of results without bias correction decreased by 14% to 23% by applying the bias correction. It can be seen that the bias correction method considering the coastline complexity can improve the model's accuracy. Finally, the effect of the proposed bias correction on projection of future SSH was examined by performing simple sensitivity tests to consider only central pressure drops to reflect the impact of climate change. By applying the bias correction, the largest potential maximum SSHs decreased by 20% from 4.52 m to 3.61 m under the FS 3. From these results, it can be seen that the proposed bias correction method can improve the model's accuracy and prevent potential overestimation of the potential maximum SSHs on the southeastern coast of the KP. The bias correction method proposed in the study can be applied to other storm surge modeling or other types of models which have results that depend on coastal topography. Figure S1. Relation between bias and distance from typhoon center to target location, and coastline complexities (left: distance, right: coastline complexities) Figure S2. Comparison of bias correction effect against