Vulnerability Analysis of Episodic Beach Erosion by Applying Storm Wave Scenarios to a Shoreline Response Model

Recently, because of the influence of climate change on sea level change, there has been growing concern regarding the erosion of beaches, which play a role in reducing the damage caused by coastal disasters. However, despite these concerns, a comprehensive understanding of the morphodynamic relationship between hazard factors and beach erosion is still lacking. Therefore, in this study, a vulnerability analysis of beach erosion was conducted by applying the shoreline response model (SLRM) of bulk model type, which identifies the physical characteristics of relevant coefficients based on the suspended sediment movement processes. To characterize wave energy incidence, storm wave scenario modeling and extreme wave analysis were conducted using wave data of 40 years on the east coast of Korea provided by the National Oceanic and Atmospheric Administration. A dimensionless mathematical function representing the storm wave scenario was proposed as a function of the peak wave height. In addition, to examine whether the beach vulnerability curve (BVC) obtained from the SLRM is valid, it was compared with the long-term shoreline observation data conducted at Maengbang Beach. For the past 9 years, sand sampling and shoreline observations were performed at Maengbang Beach about 5 times a year. However, since observations were performed in time intervals of several months, the direct comparison with model results was impossible, so a comparative analysis through statistical analysis of shoreline variability was performed. The variability of the shoreline for each reference point followed a normal distribution with a standard deviation of approximately 7.1 m. As a result of comparing the BVC results obtained from these statistical characteristics with those obtained from the model, significant similarity was shown in the high wave condition. Finally, the model was performed on two factors (mean wave height and peak wave height) which appear in SWSF and three factors (wave energy at breaking point, beach response factor and beach recovery factor) which appear in SLRM, and by analyzing the results, an approximate formula for the BVC is derived. This novel BVC approximation equation provides an intuitive understanding of the factors that affect beach vulnerability as well as their importance, and estimates the beach buffer section required to prevent coastal facilities from being damaged by erosion during a specific period. The results of this study can help limit reckless coastal development and mitigate erosion damage.


INTRODUCTION
Numerous research efforts have been made, including establishing concepts and developing models, to evaluate the vulnerability of the coastal environment in relation to coastal disasters and climate change (IPCC, 1996(IPCC, , 2001UKCIP, 2003;UNDP, 2005;UNFCCC, 2005;Han, 2006;Yook et al., 2011). Although studies often define vulnerability differently, in general, it can be characterized to be receiving a certain amount of impact from external pressure. In particular, the erosion vulnerability of coastal areas, which are prone to coastal disasters because of the effects of climate change, such as sea level changes, is emerging as a research topic of interest. However, despite these efforts, a comprehensive understanding of the hazard factors related to beach erosion and the morphodynamic processes of shoreline retreat is lacking for evaluating beach vulnerability.
Most of the understanding of episodic erosion phenomena is based on statistical analysis of beach erosion survey data or numerical model results that rely on physical parameters related to the sediment process (Mendoza and Jimeńez, 2006;Bosom and Jimeńez, 2011;Oliveira et al., 2014;Ballesteros et al., 2017;Narra et al., 2019;Anfuso et al., 2021). Data analysis has been mainly used for the purpose of finding long-term erosion rates, and short-term erosion has recently become possible with the analysis of CCTV image data. Yates et al. (2009) estimated the convergence position of the shoreline by examining changes in the beach profile as a response to changes in the incoming wave energy via long-term field measurements. In this scenario, if the external hazard is regarded as incoming wave energy and the impact factor is shoreline retreat, the beach vulnerability can be studied. However, their research was limited because the location of the shoreline that ultimately converges when a specific wave energy value was continuously applied. Recently, Park et al. (2019) suggested a methodology for analyzing beach vulnerability by analyzing shoreline data variability observed on the east coast of Korea and the extreme behavior of incident waves. However, the methodology has not been verified because of a lack of long-term wave and shoreline survey data for comparison.
To compensate for these limitations, a method that can quantitatively evaluate beach vulnerability by applying a highly reliable numerical model is required. However, the incidence of wind waves from distant seas to the shore shows an infinitely repeating phenomenon and irregular changes. Therefore, ensuring the reliability of the numerical model is important, and it is also necessary to establish a wave scenario model that considers the characteristics of wave duration. Furthermore, it is necessary to analyze how the duration of the incident wave affects shoreline fluctuations. Specifically, it is necessary to develop a methodology for extracting a vulnerability curve by analyzing long-term wave data and reflecting the influence of the wave duration to determine the maximum receding width of the shoreline. Unlike the wave input conditions used for designing port facilities, the beach response assumes that erosion proceeds according to the incidence of high waves, and when the waves calm, it recovers to its original position. However, there are few cases in which scenario models have been constructed for this purpose. Therefore, in this study, using 40-year National Oceanic and Atmospheric Administration (NOAA) data for the east coast of Korea, a scenario model formula of incident wave height and duration that is dimensionless to the maximum wave height was derived.
Several models have been developed to numerically simulate beach profile responses to identify the causes and develop preventive responses to beach erosion problems caused by storm wave incidence (Larson et al., 1999;Lesser et al., 2004;Roelvink et al., 2009;Deltares, 2018a,b). A shoreline response model of bulk model type was first proposed by Wright et al. (1985a,b). It was used to represent the rate at which the beach state changes and applied the convergence characteristics of beach conditions following storm events in situations where the mechanism of sediment movement in the breaking zone had not yet been elucidated. Later, this equation was also employed for sandbar evolution (Plant et al., 1999) and beach slope change (Madsen and Plant, 2001). Approximately 20 years later, Miller and Dean (2004) proposed a simple governing equation of the ODE type that, based on empirical evidence, relates shoreline change to shoreline position. This equation was calibrated using periodic shoreline data.
The proposed governing equations reveal that the shoreline approaches the equilibrium form for steady-state conditions at approximately exponential rates, which is consistent with laboratory investigations (Swart, 1974) and numerical results (Kriebel and Dean, 1985;Larson et al., 1999). In addition, the recent researches (Yates et al., 2009;Davidson et al., 2013;Toimil et al., 2017) established or examined a shoreline response model using long-term field observations to evaluate its practical applicability. In this study, the ODE type model proposed by Lim et al. (2021) was employed, which derived a correlation of 74% with respect to the wave input and shoreline change observation data used for the blind test of Montaño et al. (2020). Overall, two main constants are employed, the beach response factor and the beach recovery factor, and their physical properties are analyzed by evaluating the horizontal behavior of suspended sediments. However, the peak erosion width indicating the degree of beach erosion vulnerability has not been presented due to the various conditions of wave incidence yet even though the parameters a r and k r of the SLRM related to the beach property of erosion vulnerability are determined.
The main objective of this study is to obtain a relational expression for the vulnerability of the shoreline to episodic erosion according to the impact of storm wave energy given as a function of the peak wave height. Figure 1 shows the detailed research process of applying the SLRM. Section "Extraction of Storm Wave Scenario from NOAA Wave Data" introduces the process of extracting a storm wave scenario function (SWSF) by analyzing NOAA wave data for Maengbang Beach on the east coast of Korea. Section "Derivation of Beach Vulnerability Curve by the Shoreline Response Model" presents a method for deriving the vulnerability curve equation by applying this scenario function to the ODE SLRM, which has been recently proven to exist high reliability. In Section "Beach Vulnerability Curve for Maengbang Beach, " in order to analyze the feasibility of this method, the dependence of beach erosional width on wave frequency was obtained using the results of periodic shoreline surveys. In addition, the results of erosion vulnerability are presented in combination with the extreme wave analysis, which corresponds to the most fundamental vulnerability analysis method, and are compared with the results obtained from the SLRM. In Section "Discussion, " the effect of duration on erosion vulnerability was discussed by modulating the duration of the SWSF. An approximate formula for erosion vulnerability was also developed by analyzing the influence of factors on erosion vulnerability. Finally, the main conclusions are presented and future research paths are proposed based on the results of this study.

Study Site
Maengbang Beach, located in Samcheok City, Gangwondo, on the east coast of South Korea, is a linear coast stretching over approximately 4.6 km in the northwest to southeast (NE-SE) direction. The average beach width is approximately 48.0 m, and the shoreline has a mixture of straight and bow shapes. Maengbang Beach is a sandy beach near Osipcheon and Maeupcheon (Supplementary Figure 1). Supplementary Figure 2 shows its tide table and wave rose diagram (Ministry of Ocean and Fisheries, 2017). The study site has a small spring tidal range of 15.6 cm and total tidal range 33.4 cm with little tidal fluctuation, and the annual mean wave incidence at the study site has a significant wave height of 1.14 m and a significant wave period of 7.78 s. The main wave direction in the deep sea shows 39.5 • , and compared to the shoreline inclination angle of about 42 • , there is a possibility of net littoral drift to the south.

National Oceanic and Atmospheric Administration Wave Data
For the duration of wave actions, wave influence is reflected in the change of beach profiles and most storm waves show a tendency to develop, peak, and then decrease (Corbella and Stretch, 2012). In this study, wave data from NOAA was used for the storm wave scenario model analysis. The annual changes in the significant wave height and significant wave period of Maengbang Beach from 1979 to 2018 are shown in Supplementary Figures 3,4, respectively. Supplementary Figure 5 is a graph showing the annual change of annual mean significant wave height, peak period, and dominant wave direction. Supplementary  Figures 5A,B show the annual variations of values corresponding to 10% intervals from 10 to 90% from the smallest value out of about 365 × 8 (2920) data every year using wave height and period data at 3 h intervals for 40 years. And Supplementary Figure 5C shows the values corresponding to 10% intervals from 10% to 40% in the +/− direction based on 34 • N using the wave direction data, showing the change by year for 40 years. It can be seen that the annual mean wave height, period, and wave direction over the past 40 years has been maintained without significant change. Therefore, the amount of longterm change in the wave environment is negligible, and the effect of climate change in this sea area can be found to be very insignificant.
The wave height-period joint probability distribution obtained from the significant wave height-peak period data for these 40 years is shown in Supplementary Figure 6, and the corresponding spectrum is shown in Figure 2. The frequency spectrum was obtained using the joint distribution by converting the period into frequency (Sorensen, 1993) and wave height into spectrum. The directional spectrum is shown in Figure 3. Also, Figures 2, 3 show the dotted line results of the spectrum and wave direction distribution obtained at intervals of every 10 years. The peak frequency is 0.12 Hz, and the incoming deep-sea wave direction showing the peak wave energy was calculated to be 25 • clockwise from true north.

Analysis of Temporal Evolution of Wave Scenario
In order to reproduce the storm wave scenario for Maengbang Beach, the method proposed by Kim et al. (2021) is adopted in this study. Using 40-year NOAA data (38.5 • N, 129.0 • E) for Bongpo Beach, located about 111 km north of Maengbang Beach, Kim et al. (2021) normalized the wave height and period data to peak wave height values and the time scale as a function of peak wave height. Therefore, the storm scenario for Maengbang Beach was reproduced similarly using 40-year wave data of NOAA grid point (37.5 • N, 129.5 • E) closest to Maengbang Beach. Supplementary Figures 7, 8 show the wave height and wave period scenarios, respectively. The dotted lines indicate the mean and maximum values of all storm wave conditions exceeding a peak wave height in deep water (H o p ) of 4.5 m. The results show that the curves have a good fit with the extreme wave scenario of the NOAA data. In addition, the scenario function of the storm wave height affecting Maengbang Beach is expressed in terms of dimensionless time t , as given in Eq. 1. (1) Where A H = 1.28 × 10 −9 , B H = 0.00007, and t pH =−0.013. The scenario function of the storm wave period can be expressed as: (2) Where A T = 4.1 × 10 −7 , B T = 0.0042, and t pT = −0.039. The time is converted to days from the dimensionless time t according to the peak height H o p (m) as follows: (3) Supplementary Figure 9 compares the time series data wherein the peak wave height and time are dimensionless for wave scenarios with peak wave heights ranging from 4.5 to 7.5 m. These values were obtained using Eq. 1. Similarly, Supplementary Figure 10 compares the time series data wherein the wave period is dimensionless for wave scenarios with a similar range for peak wave heights. These values were obtained using Eq. 2. The dimensionless time is converted to the real time (day) via multiplication with

Wave Energy at Breaking Point
To apply the SWSF to the SLRM, it is necessary to calculate the wave height at the breaking point using the offshore wave information. If waves are induced in deep water via wave refraction and wave shoaling and they are broken near shore, the height H b at the breaking point can be obtained according to the law of conservation of energy. In a linear coastal terrain, wherein the isochore line is parallel to the shoreline, assuming that waves break in the shallow water, Eq. 1 can be rearranged in terms of H b to obtain the following: Then, by substituting h b = H b /γ, where γ is the wave breaking coefficient, the relationship can be reorganized as: where the wave breaking coefficient γ is assumed to be 0.55, which is the generally applied value for significant wave heights. The wave direction θ b at the breaking point (assumed to be in shallow water) was obtained by applying Snell's law as follows: This effect of refraction is neglected when calculating the breaking wave height, which is the input data of SLRM, from the SWSF of deep sea waves, but Eq. 6 is applied when performing the extreme value analysis on the breaking wave height in Section "Extreme Wave Height Analysis of Incidence Waves".  shows the wave scenario functions determined using the peak wave height in deep water H o p for the real-time scale t real and the corresponding wave energy at the breaking point.

Governing Equation of the Shoreline Response Model
The governing equation used to determine the maximum retreat width of the shoreline as an indicator of erosion vulnerability is given by Eq. 7. This equation was empirically obtained by Miller and Dean (2004), and recently physical validity has been verified by Lim et al. (2021) by using the horizontal behavior concept of suspended sediment. In addition, applying about 11 years of wave input data (via model prediction) and shoreline data (via CCTV analysis) provided by Montaño et al. (2020), Lim et al. (2021) achieved satisfactory model performance showing a 74% correlation.
where y is the shoreline erosion width and k r is the beach recovery factor, which is related to the recovery speed at which the suspended sediment returns to its original position after erosion.
Further, E b is the energy at the breaking point for the wave scenario model described in Subsection "Analysis of Temporal Evolution of Wave Scenario, " and a r is the beach response factor, which is a positive value in Eq. 7, and its value can be obtained from the sedimentation trend curve equation presented by Yates et al. (2009) or simply from D 50 , as suggested by Kim and Lee (2018).

Estimating Beach Response and Recovery Factors
Two physical parameters exist in the shoreline response equation (Eq. 7). The first is the beach response factor a r and the second is the beach recovery factor k r . The first factor, a r , can be easily estimated from D 50 , as described below, while the second factor, k r , can be obtained from the variation of the shoreline survey data. Based on Dean's equilibrium beach profile formula, Kim and Lee (2018) proposed a linear relationship between the equilibrium shoreline position S eq (negative value) and wave energy at the breaking point E b that is similar to the field observations of Yates et al. (2009).
where H b is the representative wave height,γ is the breaking coefficient, E b is the annual mean wave energy at the breaking point, and m i refers to the initial slope, which can be estimated as a function of beach scale factor A as follows: where A is a function of D 50 estimated from the table provided by Dean (1977). Yates et al. (2009) found that for Torrey Pines Beach, California, A has a value of f = 1.51m 1/2 for D 50 = 0.23 mm. Yates et al. (2009) previously proposed an equation similar to Eq. 8 based on Dean's equation of the equilibrium beach profile, using which they identified the shoreline erosion/accretion trend according to the wave energy and associated shoreline position of the California coast. Overall, they obtained a linear relationship between wave energy E and shoreline erosion width S eq where erosion and sedimentation are balanced: where a r is a proportional constant between the wave energy and shoreline position, and a larger value indicates lower vulnerability of the beach, and vice versa. Therefore, a r tends to increase as the grain size increases, and tends to decrease as the grain size decreases. a r is considered to be a factor indicative of protection against beach erosion because a high value mitigates beach erosion against high waves, and its reciprocal was determined to be the vulnerability proportionality constant in this study. Note that  it is assumed that b in Eq. 10 is relatively negligible for a high wave incidence. An approximate solution of the proportional constant a r (negative value) given in Eq. 8 can be obtained using the following equation: According to Eq. 11, a r = 0.0018 m when the sand particle size is 0.1 mm, and a r = 0.0111 m when the sand particle size is 1.0 mm at H b = 6.0 m andγ = 0.55 (for the significant wave). Figure 5 shows how the linear relationship between the equilibrium MSL and wave energy at the breaking point changes according to the median particle size. Note that the beach response factor a r is regarded as a positive value in Eq. 7. Figure 6 shows the wave scenario at the breaking points, and the numerical solutions of Eq. 7 for k r = 0.04, 0.06, and 0.08 d −1 . The results show that the wave input that causes shoreline change lasts for 2-3 days, but the evolution of the shoreline shows that it is affected for more than 30 days before returning to the original initial shoreline position. For k r = 0.06 d −1 , it is estimated that approximately 50 days are required for a 95% recovery after maximum shoreline erosion occurs. In this simulation, H p = 6.0 m and a r = 0.007 m were used for the peak wave height in deep water and beach response factor, respectively.
Predicting Peak Erosion Width Using the SWSF Formula Figure 7 shows the results obtained using Eq. 7 to determine the peak erosion width for peak wave heights of 4.5, 5.5, 6.5, and 7.5 m. The results of the contour line correspond to the proportionality constants of the vulnerability curves. The ranges of the two factors considered were evaluated based on whether they corresponded to a range of values corresponding to a general sand beach. The temporal change in the wave energy at the breaking point given in Eq. 7 was obtained after converting the significant wave height and peak wave period given by the scenario functions into the wave height at the breaking point.
Based on the ODE in Eq. 7, because dy/dt becomes zero at the time of peak SL retreat, when the elapsed time from the peak wave height passing to the peak erosion retreat is τ, E b (τ)/a r indicates the y peak . Therefore, y peak given as a function of τ is: Figure 8 shows the results of the τ analysis for H = 1.14 m. As shown by Eq. 13, τ is a function of only k r , regardless of a r .
τ day ∼ = −0.8 ln k r + 0.6 where H o p is the peak wave height in deep water, H o p and H have the same unit (m), and k r has a unit of d −1 . Figure 9 shows a comparison of the BVC obtained from the model results using Eq. 7 and that obtained from the approximate results using Eq. 12 for the three different values of a r and k r . Compared with each other, they exhibit satisfactory results.

BEACH VULNERABILITY CURVE FOR MAENGBANG BEACH
In Section "Derivation of Beach Vulnerability Curve by the Shoreline Response Model, " it was confirmed that the peak erosion width and peak wave energy at the breaking point had a linear correlation without substantial error. The changes in the proportional constant of erosion vulnerability based on a r and k r were calculated using the numerical model results. Therefore, if a r and k r can be determined, based on the sand properties of the beach, we can obtain the erosion vulnerability of the beach. In this section, we demonstrate the validity of this result by comparing it with the vulnerability curves obtained from analysis of the shoreline survey data observed for approximately 9 years at Maengbang Beach in Korea and the NOAA wave data near Maengbang Beach.

Shoreline Survey at Maengbang Beach
Maengbang Beach is one of the target sites for the erosion status survey of Gangwon-do. This sandy beach, with a total length of 4.6 km, stretches to the southeast, and its main incident waves originate from the northeast. Since 2010, GPS shoreline surveys have been conducted at the site four times a year. Beach profile measurements up to a depth of 15 m, sufficiently deeper than FIGURE 8 | Elapsed time τ from peak wave height to peak erosion width w.r.t the beach recovery factor k r .
the closure depth, were obtained twice a year at 150 m intervals, and the median particle size survey was conducted in the swash zone. The beach profile was measured by potable GNSS on foot from the survey reference point to a water depth of 1.5 m and from 1.5 m to 15 m water depth was surveyed using a bathymetry survey boat equipped with an echo-sounder (single beam). The GNSS model used for shoreline surveying is GX1230 with an accuracy of H: 3 mm + 0.5 ppm. The echo-sounder used for the bathymetric survey is a single beam, the model name is AquaRoller 200S, and the measurement range is 5 ∼ 200 m. These datasets were utilized as basic data to evaluate the coastal erosion grades (A, B, C, and D). In the event of the erosion of the coast, the datasets serve as basic data for cause analysis. Figure 10 shows the 13 reference points for the survey conducted at Maengbang Beach, which was divided into four survey zones, wherein the mean shoreline and erosion control line (erosion limit for 30 years) were obtained according to the statistical characteristics of the shoreline data. These two shorelines were used for (1) evaluating the management status of Korean coastal regions in terms of conservation and disaster prevention, and (2) evaluating whether the layout design of soft and hard engineering structures of the coastal improvement projects implemented to reduce coastal erosion are achieving the management target lines (mean target shoreline and erosion prevention line).
The average D 50 obtained for each zone and the corresponding beach scale factor and beach response factor are listed in Supplementary Table 1. Although the average sand particle diameter of Zone 4 was relatively small, the average particle diameter of the entire target beach was 0.656 mm. On the eastern coast of Korea, the net littoral drift flows from the north to the south according to the analysis of long-term wave data and the angle of the coastline (Kim et al., 2001). Owing to this littoral drift environment, the sand on the southern beach may have a finer grain size, but the difference is expected to be small.

Variability Analysis Using Shoreline Survey Data
The shoreline survey of Maengbang Beach, the subject of this study's beach response analysis, was conducted about five times a year since 2010 as a major observation item in the coastal erosion survey, and was measured at 13 reference points. In the erosion status survey, the shoreline survey is, in principle, conducted once a season. However, in some cases, shoreline surveys were more carried out once or twice a year. Therefore, it resulted in an average of about five shoreline surveys. Figure 11 shows the distribution histograms for each survey zone using 192 shoreline survey datasets collected at 13 reference points (see Figure 10),  for which measurements were conducted 47 times from 2010 to 2018. The figure illustrates the probability distribution of the observed shoreline data with reference to the mean values for each zonal reference point. These results were compared to a Gaussian distribution and were found to exhibit high similarity. The statistical frequency analysis of these shoreline observations allows researchers to assess the risk of beach erosion.
Zone 1 is located on the northern side of Maengbang Beach, and Zone 4 is located on the southern side of the beach. The results of each zone reasonably follow a normal Gaussian distribution, and according to their statistical characteristics, the monitored erosion width (MEW) can be obtained for each return period (yr). The standard deviation σ for each zone is listed in Supplementary Table 2, as well as the MEW for each zone, which was calculated from the mean value for each return period based on the frequency analysis. These results represent erosion damage by frequency and can therefore be referred to as a risk curve. Conversely, if the extreme wave height analysis result of the NOAA wave data is applied instead of the return frequency, the coastal erosion vulnerability curve (damage curve according to hazards) can be derived.
To verify the normality of the observed shoreline variations, a chi-squared goodness of fit test was performed, confirming that the results follow a normal distribution at the 1% significance level. A qq plot is shown in Supplementary Figure 12. Supplementary Table 2 lists the standard deviation σ of the shoreline survey data for each zone and the MEW for each return period, which was obtained from a variability analysis of the data. The larger the standard deviation, the greater the variability, and thus, the greater the erosion width. Zone 4, which consists of fine sediment, exhibited the highest MEW.

Extreme Wave Height Analysis of Incident Waves
For the extreme wave height analysis, the 40 years long-term wave estimation data provided by NOAA, consisting of 3 h intervals from January 1979 to December 2018, were used. The NOAA deep-sea wave data were converted by applying the shoaling, refraction, and significant wave breaking conditions of Eqs. 4-6 to obtain the wave height at the breaking point. Supplementary Figure 13 shows the results of the extreme wave height analysis conducted using the Gumbel method, and Supplementary Table 3 summarizes the extreme wave height at breaking point H F obtained for each frequency F, which is the reciprocal of the return period. Approximately half of the wave direction components that could not approach the target beach because the wave direction was too wide were excluded from the extreme distribution analysis. The relationship between the extreme wave height (m) and frequency F (yr 1 ) is as follows:

Analysis of Erosion Vulnerability for Maengbang Beach
As mentioned in Section "Variability Analysis Using Shoreline Survey Data, " the MEW for each return period was obtained via an analysis of fluctuations in the monitored shoreline. Therefore, using Eq. 14 to obtain the extreme wave height for each return period, an beach vulnerability curve showing the erosion width for each wave height could be obtained. Supplementary Table 4 lists the calculated MEW for each zone of Maengbang Beach for each extreme wave height at a breaking point of 4.5 m or higher. Figure 11 compares the BVC directly obtained from the extreme breaking wave height and the MEW listed in Supplementary  Table 4 with the BVC numerically obtained by applying the SWSF to the SLRM. In the SLRM, the beach response factor a r was obtained using Eq. 11 using the sand particle size data, and the beach recovery factor k r was applied as the value with the best fit, as compared to the MEW data for E b > 2.72 m 2 . The k r values obtained for zones 1-4 were 0.0698, 0.0624, 0.0552, and 0.0796 d −1 , respectively. When H b > 6.6 m, at which E b > 2.72 m 2 , there is a satisfactory trend, but when H b < 6.6 m, the values interpreted from the observations are lower than the model predictions, and unlike the model predictions, they do not converge to zero. Instead, the values exhibit a convergence to E b = 1.73 m 2 , corresponding to the 1-year return period. When the SWSF is shorter than the 1year return period, the selection of a curve along the maximum coverage line as the scenario model may lead to discrepancies, as shown in Supplementary Figures 9A, 10A. Zones 1 and 4, which are located at the ends of the beach, show slightly larger k r values than zones 2 and 3. It is presumed that fitting with a larger k r value is probably because of the greater variability at the beach ends, mainly because of the influence of seasonal littoral sedimentation. While k r can be determined from the particle size of the sand, estimating its value is difficult and requires additional research with on-site monitoring. Note that this study is limited to analyzing beach vulnerability in the state where a r and k r are estimated.

Comparison With Beach Vulnerability Curve From Shoreline Response Model
As shown in Figure 11, zones 1 and 4, unlike the zones located in the center of the beach, had relatively high k r values compared to the sand grain size, which may be related to the occurrence of littoral drift due to the incidence of oblique waves at the beach ends. On this basis, considering the loss of suspended sediment due to littoral drift, Eq. 7 is modified as follows: where ε corresponds to the loss rate of suspended sediment caused by the occurrence of littoral drift. If ε is greater than 1, it is considered that suspended particles are introduced from the outside, and when ε is less than 1, it is considered that particles are discharged to the outside. Figure 12 shows how the erosion width varies with ε. In particular, when ε = 0.2 is applied to zones 1 and 4, a reasonable k r is obtained, as compared with that in the surrounding zones. The k r values obtained for zones 1-4 were 0.0628, 0.0624, 0.0522, and 0.0698 d −1 , respectively.

DISCUSSION
In this study, the SWSF and SLRM were applied to determine the BVC using only cross-shore beach sedimentation based on the values of the beach response factor a r and beach recovery factor k r . The BVC may greatly vary depending on the duration of the coast-specific input SWSF, as well as a r and k r . Specifically, even with the same peak wave height, the erosion width may vary according to the SWSF when subjected to a limited duration. Therefore, in this section, we attempt to determine the time constraint characteristics of the BVC that is less than equilibrium owing to the effect of a limited duration, as compared with the extreme BVC, which has an equilibrium erosion width and the duration of the incoming wave is infinite. In general, we attempt to determine the time constraint characteristics of a BVC that is not in equilibrium by comparing it with the extreme BVC. The extreme BVC corresponds to the results obtained from the field experiments of Yates et al. (2009). However, note that the resulting erosion width is too high to be realistic, and the actual occurrence of erosion of this magnitude would be catastrophic.
When the equilibrium state is reached, the beach erosion vulnerability is roughly expressed by the wave energy at breaking point E b and the beach erosion vulnerability factor f v : where ESLP is the equilibrium shoreline position and f v is the reciprocal of the beach response factor a r . This equation is based on Eq. 8, which was modified by ignoring E b , as it was considered to be relatively small under high wave conditions. Overall, as the vulnerability factor increases, more beach erosion occurs, indicating that the beach erosion vulnerability increases. However, the ESLP refers to a virtual erosion state formed by consistant waves incoming with an infinite duration. Unlike the ESLP, the duration-limited shoreline position (DSLP), which is affected by duration and has a smaller erosion width, can be expressed as follows: Therefore, the ratio of the erosion width actually generated by the SWSF compared to the equilibrium erosion width can be defined as the shoreline position ratio µ, as follows: Supplementary Table 5 lists the calculated τ according to k r and H o p , revealing an elapsed time of approximately 1.5-3 days. Supplementary Table 6 lists the changes in the µ value defined by Eq. 18. These results are valid for the wave conditions at the Maengbang coast. Unlike the ESLP, which is the result of infinite waves of a certain height, the erosion width of the beach that occurs under the actual duration-limited condition, the DSLP, shows a much smaller erosion width, exhibiting values of only approximately µ = 0.02−0.11. This is because of the characteristics of the real-world incident storm wave scenario, wherein wave energy is not continuously applied at a constant value but increases and then decreases, not allowing an equilibrium state to be reached. Figure 13 shows H has the effect of spreading SWSF laterally. Overall, as H increased, the influence of the duration increased. Although H, which is defined as the mean wave height, is unlikely to be greater than 2 m, we assessed the effect of large H values to observe the effect of duration. Figure 14A shows that µ becomes larger owing to the effect of longer duration, at constant values of a r and k r . Because µ tends to decrease as wave height increases, H o p /H was applied to eliminate this change and obtain a result that barely changes not only with the wave height but also the H, as shown in Figure 14B. In addition, considering the change in k r , we obtained the following equation: where a r and k r have units of m and d −1 , respectively, and the three coefficients α 1 , α 2 , and α 3 have values of 1.5, −0.275, and 0.911, respectively. By inserting Eq. 19 into Eq. 17, we obtain the following equation, which is useful for the practical application of the DSLP: Thus, the peak erosion width can be conveniently estimated by using Eq. 20 if the incident peak wave information and characteristic coefficients of the beach responses a r and k r are determined. This gives an idea of how wide the buffer zone should be to avoid the damage to the hinterland for the return period of interest. However, this result ignores the effect of longshore sediment transport and only considers the cross-shore sediment transport process. We examined whether Eq. 20 provides a satisfactory solution, as compared with the SLRM result. Figure 15 shows the results of the comparison under the same conditions as in Figure 9, revealing fairly agreeable results.
The reliability of the BVC (Eq. 20) obtained from the shoreline response model was reviewed by comparing the results obtained  from the SLRM by applying SWSF with the results obtained from the statistical analysis of periodic shoreline survey data. All of the 47 shoreline survey data taken in the target area are shown in Supplementary Figure 14, and the average by statistical analysis and the erosion width by return period are plotted thereon. The applicability of the BVC proposed by Eq. 20 was verified and it was confirmed that it showed satisfactory reliability. Supplementary Figure 15 shows the results obtained from BVC among the results shown on the three dimensional LiDAR image. The beach recovery factor k r is the result of considering the effects of littoral drift as described in Section "Comparision with BVC from SLRM".

CONCLUSION
Vulnerability is measured as the degree of damage to hazard intensity. In this study, therefore, the storm wave incident scenario to the target beach expressed in terms of the peak wave height, which is the hazard intensity, was established and the scenarios were applied to the SLRM, and the vulnerability curve was extracted by correlating the result of the erosion width, which is the degree of damage. This methodology can be applied to beaches where there are few tidal ranges. Although it may be slightly affected by other factors, such as berm height and initial slope, these effects are considered to be insignificant. The storm wave scenario function (SWSF) obtained by analyzing NOAA wave data was used as the input data for the SLRM. The numerical results of the model provided satisfactory results compared with the results of shoreline observations, which were conducted five times a year for 9 years on the eastern coast of Korea. Further, while the model did not provide satisfactory results for wave heights with a return period of approximately 1 year or less, it showed fairly good agreement under high wave conditions, which is more meaningful for vulnerability analyzes.
By analyzing the factors affecting erosion vulnerability via the SLRM, an approximate equation with very good consistency was developed. This equation is given by the peak wave height at deep water H o p and the mean wave height H in relation to the SWSF, as well as the beach response factor a r and beach recovery factor k r in relation to the SLRM. However, although it is not very difficult, the H o p needs to be converted to H b by considering wave shoaling, wave refraction, and wave breaking. This BVC approximation provides an intuitive understanding of the factors that influence beach vulnerability and estimates the length of the beach buffer zone required to prevent erosion damage to hinterland facilities over a specific return period. This is expected to provide essential information for limiting reckless coastal development and mitigating erosion damage.
Herein, the SWSF was determined by analyzing the storm wave scenarios with wave heights of 4.5 m or higher from NOAA wave data, which included long-term wave estimation data over 40 years, at 3 h intervals from January 1979 to December at 38.5 • N, 129.0 • E near Maengbang Beach. To examine the validity of the BVC obtained from the model results, shoreline surveys were conducted approximately five times a year for 9 years and sand size data were used. In the range of E b > 2.72 m 2 , the results of SLRM and statistical analysis of MEW variability show quite similar patterns, so the prediction of the peak erosion width by SLRM combined with SWSF is considered reasonable. The beach recovery factor k r in SLRM was applied as the best fit value compared to the MEW data. The k r value is a physical factor affecting the recovery of the shoreline to its original state, and this value is expected to be related to the grain size of the sand. Future research on this is needed.
Because the tidal wave difference is small along the eastern coast of Korea, the SLRM can be directly applied to the coast. However, in the case of the western coast of Korea and other coastal areas with high tidal variations, the influence of incoming high waves evenly affects the intertidal zone, making it difficult to directly apply the methodology of this study. Therefore, in the future, the SLRM must be improved to assess the erosion vulnerability in coastal areas with high tidal variations, to reflect the effect of tides.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
T-KK: conceptualization, acquisition of data, visualization, and writing -editing. CL: acquisition of data, visualization, and writing -editing. J-LL: conceptualization, supervision, validation, and writing -original draft and review. All authors contributed to the article and approved the submitted version.