Developments of Dynamic Shoreline Planform of Crenulate-Shaped Bay by a Novel Evolution Formulation

: In this paper, a simple evolution formulation, based on polar coordinate, is proposed to efﬁciently and accurately simulate the dynamic movement of sandy shoreline in a crenulate-shaped bay. The consistent actions by seasonal waves and swell usually result in severe erosion or deposition along sandy shoreline, so some mathematical formulations and numerical models, based on the concept of headland control, have been recently proposed to forecast and protect the shoreline planform. A simple and general formulation is derived in this study by considering the balance of longshore sediment transport, since the accurate prediction of dynamic movements of high planform-curvature shoreline in the shadow zone behind a headland is essential and critical to the headland control. The proposed formulation can be directly adopted to accurately and simply simulate the temporal variations of shoreline in both of the hooked zone, where is protected by a headland, and the unhooked zone, where is straight attacked by incident waves. Numerical results and comparisons of temporal variation of shoreline between two headlands are provided in this paper to demonstrate the accuracy and efﬁciency of the proposed evolution formulation. Besides, different time increments and different numbers of control volume are adopted to examine the merits of the proposed numerical model.


Introduction
At part of the junction of land and sea, there are some beautiful sandy beaches, which can naturally provide leisure spots or coastal protections.Since the shoreline of sandy beaches are formed by the combinations of complicated natural forces, such as wave, current, tide, sand, wind, etc., it can be easily observed that the planform of crenulateshaped bays will dynamically deform due to shoreline continuous retreat and advance from these natural forces.In order to protect the coastal zone and remain the leisure spots, some methods of coastal protection have been proposed in the past decades.
Jetty, groin and breakwater are the most common structures in coastal engineering.These structures can maintain the shoreline and present the effect of crenulate-shaped bays.Özölçer [1] experimentally studied the relationship between jetty and shoreline in laboratory and received some innovative results.Putro [2] used parabolic bay shape equation (PBSE) to analyze the effectiveness of groins in practical application at Nusa Dua Beach.Although these jetty and groin can effectively maintain the beach, the up-coast and down-coast imbalance caused by jetting effect may lead to the erosion of the originally safe coast due to the lack of sand sources.Liang [3] studied the phenomenon of coastal erosion down-coast once the groins on the coast of Yilan at Taiwan was built.From the viewpoint of coastal protection, although these jetty and groin can effectively trap the longshore sediment transport, they cannot provide sufficient protection for the coastline and cannot effectively reduce wave energy in the face of typhoons.
In 1952, Penney and Price [4] proposed the diffraction theory of waves by considering a semi-infinite breakwater.The shore protection manual [5] drew the distribution of wave height once the incident waves were diffracted by the breakwater.After that, many scholars successively proposed different types of breakwaters and proved that these breakwaters can provide effective coastal protection [6,7].In addition, the static and stable waters behind the breakwater cause the drifting sand to settle and begin to silt up, so many scholars have begun to use this engineering method to maintain the beach.For example, Laustrup [8,9] and Zhu [8,9] investigated the relationship between offshore breakwater and salient.While protecting the coastline, the beach behind the offshore breakwater can meet the needs of modern people for hydrophilic space and landscape, so the theory of artificial headland has been developed [10].
In the headland control, an artificial structure or a natural rocky island is placed at sea and close to sandy beaches.Due to the continuous movements of sediment by natural forces, a beautiful sandy planform can be established between two headlands or a headland and a control point, such as tombolo.Although the sandy beaches of crenulate-shapes bays can be formed by using the headland, the unsteady weather and wave climate will randomly move the sediments, which may result in unsteady erosion or deposition along shoreline.Thus, beach nourishment is sometimes necessary to protect the shoreline in advance.From the above discussions, it is noticed that the shape and the movement of shoreline in sandy bays play an important role in coastal protection.As a result, some methods, which can simulate the shape of shoreline and its possible dynamic movements, are essential for practical engineering design.In this paper, we proposed an evolution formulation, based on polar coordinate, to simply and accurately predict the temporal variations of shoreline planform of a crenulate-shaped bay.
The shape and movements of shoreline of sandy beaches is the integrated results by wave climate, current, wind, etc., so in the past half century many researchers paid much attention on the correct descriptions of the steady sandy planform of embayed beach.It is found by [11] that the planform of embayed beach in static equilibrium is mainly dependent on the obliquity of incident wave and the control line.The shape of shoreline can be described by a parabolic bay shape model, which leads to various study of the headland control.Later, Gonzalez and Medina [12] modified the parabolic bay shape model, proposed by [11], by considering the analytical model of null drift current, so their proposed model can be adopted to predict the planforms of undeveloped manmade and natural beaches.Their proposed methodology has been successfully applied to describe the shapes of sandy beaches along Atlantic and Mediterranean coasts of Spain.In addition, Iglesias et al. [13] adopted the artificial neural networks (ANN) to study the planform of embayed beaches in static equilibrium.Three ANN models, which are all based on the feed-forward backpropagation networks, have been proposed to study the shape of crenulate-shaped beaches.Since the training data and the testing data are acquired from the planform of real beaches, their results are compared very well with some conventional planform models.Hsu et al. [14] provided a comprehensive review for the past study of static bay beach, since the scientists and engineers of costal engineering have studied the steady planform of embayed beaches for more than half century.The developments of empirical models and the methodology for stability verification of a headland-bay beach are clearly discussed in the review.The applications of the relevant study for coastal management are also provided.Besides, according to plenty researches of the PBSE, some computer software have been established and presented.Among them, the MEPBAY and the SMC software are the most famous and widely used for engineering design.Raabe et al. [15] completely compared two of them from the theoretical model to the accuracy of numerical results.Considering the required accuracy of solutions and the computational efficiency, it is suggested in their study that the MEPBAY can be used in the pre-design phase, while the SMC is more suitable for the design phase.Although Water 2022, 14, 3504 3 of 22 some useful formulations [11,12,16,17] and computer packages have been proposed and validated to predict the static equilibrium of shoreline planform, it is worthy to develop a model to accurately simulate the dynamic movements of shoreline planform due to continuous natural forces.
In addition to the developments of static-equilibrium planform, some study in the past decades have been presented to accurately simulate the temporal change of bay shape.Pelnard-Considere [18] proposed a one-line formulation for shoreline change by considering the retreat and advance, induced by the longshore gradient in longshore sediment transport.Following the one-line formulation proposed by [16,18,19] described the well-known numerical model, called the GENESIS, to simulate the dynamic change of shoreline.The GENESIS, which includes wave shoaling, refraction and diffraction, can calculate the shoreline change due to groins, jetties, detached breakwater, seawall and beach fills.But, the governing equation in the GENESIS is discretized by using the Cartesian coordinate, which may pose some limitations to the GENESIS.In order to improve the mathematical formulation in the GENESIS, Weesakul et al., (2010) [20] derived the one-line formulation by considering the polar coordinate in the hooked zone, where is in the lee of a headland, and the Cartesian coordinate in other area.Moreover, Kaergaard and Fredsoe (2013) [21] proposed another novel one-line formulation to model the dynamic evolution of shorelines with large curvature.Their proposed formulation is derived by using local orthogonal coordinates, which is parallel and normal to every shoreline segments.The accuracy of numerical results in their model is strictly depended on the spectral wave field, the longshore current flow field and the sediment transport field, which should be acquired by using a commercial software, the Mike21FM.Though some successful methods have been proposed to predict the dynamic movement of shoreline planform, in this paper we proposed an evolution formulation to simply and accurately simulate the transition state of shoreline planform in a crenulate-shaped bay [16,[18][19][20][22][23][24][25].
In this paper, we used only the polar coordinate to develop the evolution formulation, which is extended from the study by [16,20].Weesakul et al. [20] developed a numerical model to simulate the transition state of shoreline planform of a crenulate-shaped bay by utilizing both of the Cartesian and the polar coordinates.The polar coordinate is adopted along the shoreline in the shadow zone of the headland, while the Cartesian coordinate is used along the shoreline in the unhooked zone, where is not protected by the headland.Although the numerical results are in good agreements with other experimental and numerical solutions, there are some parts of their models, which are deserved to be improved.For example, the treatment along the interface between two different coordinates may possibly lead to instability of numerical simulation.Besides, to calculate the quadratic term in the polar coordinate-based evolution formulation is quite expensive in the viewpoint of computational efficiency.In order to simplify and improve the numerical model, proposed by [20], we utilized only the polar coordinate to develop an evolution formulation by considering the mass conservation of sediment transport and a rectangular shape of control volume.From the above discussions, it is apparent that the evolutionary formulations, proposed by the GENESIS [16], Weesakul et al. [20] and the present paper, adopted different coordinate systems to describe the shoreline.As a result, the directions of shoreline advance and shoreline retreat are very different from each other.In the GENESIS and in the unhooked zone in [20], the direction of shoreline advance and shoreline retreat are identical to the y direction, which is perpendicular to the x direction in Cartesian coordinate.On the other hand, in the hooked zone in [20] and in the present paper the direction of shoreline advance and shoreline retreat is equal to the radial direction in polar coordinate.
The dynamic evolution of shoreline planform between two headlands are considered in this paper.By placing the origin of the polar-coordinate system at the tip of one of these two headlands, we can describe the spatial positions of sandy shoreline from one headland to the other headland by using only the polar coordinate.The considered shoreline can be divided into several control volumes, so we can derive the evolution formulation by using a plan rectangular control volume and the conservation law of mass.Since the quadratic term in the evolution formulation is derived by using a trapezoidal control volume [20], in this paper we adopted the rectangular control volume to avoid the quadratic term and simplify the one-line formulation.For the derivation of the one-line evolution formulation in Section 2, only the longshore sediment transport is considered.It is assumed that there is no longshore sediment transport at two headlands.After the effectiveness of the evolution formulation is validated in Section 3, we first consider the change of the shoreline under the condition of volume conservation, and then add the consideration of the cross-offshore drift factor in Section 4. Accordingly, a very simple evolution formulation to accurately describe the dynamic movements of sandy shoreline is derived in this study.In this paper, the numerical results are compared with other numerical solutions to demonstrate the accuracy and applicability of the proposed formulation.

Evolution Formulation for Dynamic Shoreline Planform
One of the most famous one-line shoreline models bas been published by GENESIS, and its basic assumptions are listed as follows: shape of beach profile remains constant; shoreward and seaward depth limits of the profile are constant; sand is transported alongshore by the action of breaking waves; the detailed structure of the nearshore circulation is ignored; there is a long-term trend in shoreline evolution.
In this section, the proposed novel evolution formulation for dynamic shoreline planform is derived by using the conservation law of mass, which means does not consider line source or sink of sand in this section, and the longshore sediment transport.The proposed formulation used the polar coordinate in both of the hooked zone and the unhooked zone, so the proposed formulation can be used along the entire shoreline.The schematic diagram of the shoreline planform and two headlands are presented in Figure 1.In addition, the definition of the polar coordinate is demonstrated in Figure 1 and the origin of the coordinate system is placed in the tip of the upcoast headland.Since the origins of the polar coordinate and the Cartesian coordinate are employed at the same control point, these two coordinate systems can be represented by the following equation, where (x, y) and (r, θ) are the space position of the Cartesian coordinate and the polar coordinate, respectively.In Figure 1, it can be found that the shoreline from the upcoast headland to the down-coast headland are divided into several segments.{r i } m i=1 is the distance between the central node of the {i} th shoreline segment and the tip of the control point of the upcoast headland, as well as m is the number of total nodes along the entire shoreline.
The schematic diagram for the volume change of shoreline segment between two time steps is depicted in Figure 2.While only the longshore sediment transport within a time increment, ∆t = (t n+1 − t n ), is considered, the spatial position of the i th shoreline segment at the n th time step, r n i , is updated to the spatial position of the i th shoreline segment at the (n + 1) th time step, r n+1 i .∆A means the changed plan area by shoreline retreats or advances and D s is the closure depth.In the proposed evolution formulation, the closure depth D s is set as a constant value, which is suggested by [26].One of the basic assumptions of the one-line model is that the beach profile is constant.In addition, Q n i is the longshore sediment transport between the i th shoreline segment and (i − 1) th shoreline segment, as well as Q n i+1 is the longshore sediment transport between the i th shoreline segment and (i + 1) th shoreline segment which is along the shoreline.The directions of the two vectors of longshore sediment transport in Figure 2 are the positive direction of longshore sediment transport; meanwhile, the difference between these two arrows within one time step will result in the shoreline advance or the shoreline retreats.Besides, these two arrows should be perpendicular to the edge faces.The four boundary segments enclosing ∆A formed the control volume of the i th shoreline segment between the n th time step and the (n + 1) th time step.From Figure 2, it is obvious that the volume change of the i th shoreline segment between two time step, ∆V, can be demonstrated as follows, ∆V = ∆Q∆t (2) Wa ter 2022, 14, x FOR PEER REVIEW 5 of 22 The schematic diagram for the volume change of shoreline segment between two time steps is depicted in Figure 2.While only the longshore sediment transport within a time increment, ∆ = (  +1 −   ), is considered, the spatial position of the i th shoreline segment at the  ℎ time step,    , is updated to the spatial position of the i th shoreline segment at the ( + 1) ℎ time step,    +1 . ∆A means the changed plan area by shoreline retreats or advances and   is the closure depth.In the proposed evolution formulation, the closure depth   is set as a constant value, which is suggested by [26].One of the basic assumptions of the one-line model is that the beach profile is constant.In addition,    is the longshore sediment transport between the  ℎ shoreline segment and ( − 1) ℎ shoreline segment, as well as  +1  is the longshore sediment transport between the  ℎ shoreline segment and ( + 1) ℎ shoreline segment which is along the shoreline.The directions of the two vectors of longshore sediment transport in Figure 2 are the positive direction of longshore sediment transport; meanwhile, the difference between these two arrows within one time step will result in the shoreline advance or the shoreline retreats.Besides, these two arrows should be perpendicular to the edge faces.The four boundary segments enclosing ∆A formed the control volume of the  ℎ shoreline segment between the  ℎ time step and the ( + 1) ℎ time step.From Figure 2, it is obvious that the volume change of the  ℎ shoreline segment between two time step, ∆V, can be demonstrated as follows,  The GENESIS used the Cartesian coordinate and a rectangle to approximate ∆A, so a very simple one-line evolution formulation is derived.However, the Cartesian coordinate-based evolution formulation is not useful for high planform-curvature shoreline.In order to tackle these problems, Weesakul et al. [20] proposed a hybrid formulation.In the The GENESIS used the Cartesian coordinate and a rectangle to approximate ∆A, so a very simple one-line evolution formulation is derived.However, the Cartesian coordinatebased evolution formulation is not useful for high planform-curvature shoreline.In order to tackle these problems, Weesakul et al. [20] proposed a hybrid formulation.In the hooked zone, the polar coordinate is adopted and a trapezoid is used to approximate ∆A.On the other hand, in the unhooked zone, the Cartesian coordinate-based evolution formulation, proposed by GENESIS, is adopted.In the hybrid formulation [20], the interface treatment between two coordinates may result in numerical instability.Furthermore, the quadratic form in the polar coordinate-based formulation in the hooked zone may require more computational resources for simulation.In order to remain the simplicity of evolution formulation, and the applications of evolution formulation to high planform-curvature shoreline, a novel evolution formulation for dynamic shoreline planform of crenulateshaped bay is derived in this paper by using the polar coordinate and a rectangle to approximate ∆A.Thus, the control volume in Figure 2 is approximated by a rectangle and ∆A can be expressed as follows, In Equation ( 4), we used rectangle formula to calculate area instead of the sector control volume in [20].The first term in the right-hand side of Equation ( 4) represented the approximation of the length of an arc, which represents the length of the rectangle.On the other hand, the second term (r n+1 − r n ) of Equation ( 4) represents the width of the rectangle.The approximations of the length and the width of the rectangular control volume are depicted in Figure 3.
Wa ter 2022, 14, x FOR PEER REVIEW 7 of 22 The novel evolution formulation proposed in this paper can be derived from Equation ( 5) and can be shown as follows, Equation ( 7) is the proposed one-line evolution formulation for dynamic shoreline planform of crenulate-shaped bay.Once the physical conditions at the  ℎ time step are known, the spatial position of shoreline at the ( + 1) ℎ time step can be simply calculated by using Equation (7).Suitable initial condition for shoreline should be provided in the beginning of numerical simulation.The temporal marching scheme is identical to the explicit Euler method, since the space position of shoreline at the ( + 1) ℎ time step can be directly obtained from the physical quantities at the  ℎ time step.Thus, it can be deduced that we should use smaller ∆ so as to ensure the stability of the numerical simulation.Furthermore, the proposed formulation can be used to simultaneously simulate the change of shoreline at both of hooked zone and unhooked zone.In addition to the initial condition, two boundary conditions should be given at the two ends of shoreline.No longshore sediment transport,  1 =  +1 = 0 , is assumed, since two headlands are placed at the two ends of the shoreline.
It is evident that, from the viewpoint of geometry, the plane shape of control volume within the hooked zone is similar to a sector of annulus or a trapezoid as well as the plane shape of control volume within the unhooked zone is a parallelogram if the polar coordinate is adopted to describe the entire shoreline.In this paper, we used the rectangular control volume, which is presented in Figure 3, to derive the evolution formulation for dynamic shoreline planform instead of the trapezoid and the parallelogram.It could be expected that the inexact descriptions for shape of control volumes may possibly induce some numerical errors.However, by carefully examined the control volume in Figure 3, The rectangular control volume, which is used in this paper and described by using the polar coordinate, is depicted in Figure 3. Substituting Equation (4) to Equation (3) can lead to the following equation, The novel evolution formulation proposed in this paper can be derived from Equation ( 5) and can be shown as follows, Equation ( 7) is the proposed one-line evolution formulation for dynamic shoreline planform of crenulate-shaped bay.Once the physical conditions at the n th time step are known, the spatial position of shoreline at the (n + 1) th time step can be simply calculated by using Equation (7).Suitable initial condition for shoreline should be provided in the beginning of numerical simulation.The temporal marching scheme is identical to the explicit Euler method, since the space position of shoreline at the (n + 1) th time step can be directly obtained from the physical quantities at the n th time step.Thus, it can be deduced that we should use smaller ∆t so as to ensure the stability of the numerical simulation.Furthermore, the proposed formulation can be used to simultaneously simulate the change of shoreline at both of hooked zone and unhooked zone.In addition to the initial condition, two boundary conditions should be given at the two ends of shoreline.No longshore sediment transport, Q 1 = Q m+1 = 0, is assumed, since two headlands are placed at the two ends of the shoreline.
It is evident that, from the viewpoint of geometry, the plane shape of control volume within the hooked zone is similar to a sector of annulus or a trapezoid as well as the plane shape of control volume within the unhooked zone is a parallelogram if the polar coordinate is adopted to describe the entire shoreline.In this paper, we used the rectangular control volume, which is presented in Figure 3, to derive the evolution formulation for dynamic shoreline planform instead of the trapezoid and the parallelogram.It could be expected that the inexact descriptions for shape of control volumes may possibly induce some numerical errors.However, by carefully examined the control volume in Figure 3, it should be emphasized that the length of the control volume, (r n i ∆θ), is in the order of tens to hundreds meters while the width of the control volume, (r n+1 − r n ), is in the order of 0.1 to 0.01 m for practical engineering applications.The huge difference between the length and the width of the control volume forms a high-aspect-ratio rectangular control volume.Since the control volume in both of hooked zone and unhooked zone is extremely slender, the numerical errors, which are resulted from the inexact description for shape of control volume, should be unapparent.Therefore, in this paper we adopted the rectangular control volume and the polar coordinate in order to simply derive the evolution formulation for dynamic shoreline planform.
It can be observed in Equation ( 7) that the movements of shoreline position are mainly determined by considering the change of longshore sediment transport.In one-line shoreline model, longshore sediment transported is driven by the action of breaking waves.In this paper, we used the following formulation to describe the longshore sediment transport, which was proposed and improved by [27] and [16,28] as well as is demonstrated in Equation ( 8).K 1 , K 2 are longshore sand transport calibration coefficients, and 0.5 K 1 < K 2 < 1.5K 1 is recommend by [28].In this study, These two coefficients K 1 , K 2 are used as the same as those in [20] where tanβ means the seabed average slope, ρ s is the density of sand, ρ is the density of sea water and λ is the porosity of sand.H b is the wave height at breaking-wave point.C g = gh b is the group wave velocity and h b means the depth at breaking-wave point.g is the gravity.α bs is wave angle at breaking wave point to shoreline.The schematic diagram for α b , α s and α bs are defined in Figure 4. Substituting Equation (8) to Equation (7) will form the proposed evolution formulation for dynamic shoreline planform.The first term in the right-hand side of Equation ( 8) means the longshore transport rate owing to the obliquely incident waves; meanwhile, the second term in the right-hand side of Equation ( 8) expresses the longshore sand transport rate induced by the longshore variation of breaking-wave height.
will form the proposed evolution formulation for dynamic shoreline planform.The term in the right-hand side of Equation ( 8) means the longshore transport rate ow the obliquely incident waves; meanwhile, the second term in the right-hand side of E tion (8) expresses the longshore sand transport rate induced by the longshore variat breaking-wave height.

Verification for the Consistency and the Stability of the Proposed Model
In order to examine the consistency and the stability of the proposed model, se numerical experiments are implemented in this section.The schematic diagram of th merical example is described in Figure 5. Two breakwaters are placed in front of the i shoreline and the distance between these two breakwaters is  0 = 120.The length of

Verification for the Consistency and the Stability of the Proposed Model
In order to examine the consistency and the stability of the proposed model, several numerical experiments are implemented in this section.The schematic diagram of the numerical example is described in Figure 5. Two breakwaters are placed in front of the initial shoreline and the distance between these two breakwaters is G 0 = 120.The length of each breakwater is B = 60, while the distance between the breakwaters and the parallel initial shoreline is S = 50.From Figure 5, it is obvious that the initial shoreline can be described by three straight lines.The lengths of these three straight lines are 50 m, 180 m and 50 m, respectively.In addition, the origin of the polar coordinate is placed at the tip of the left breakwater.When the initial shorelines and some initial conditions are provided, the dynamic shoreline planform at different time steps can be simulated by using the proposed evolution formulation, Equation (7).From numerical experiment, it is found that the change of shoreline position is negligible when time is large than 50 years.Therefore, the numerical results at t = 50 years will be recognized as the steady-state solutions.In this section, the steady numerical results by using different numbers of control volume, m, and different time increments, ∆t (s), are demonstrated and compared with each other so as to validate the consistency and stability of the proposed model.
An ideal numerical example is considered and the initial shoreline is depicted in Figure 5.The initial shoreline is constructed by three straight lines.Therefore, there will be two corners, which are the discontinuities of the geometry.While the novel evolution formulation for dynamic shoreline planform is derived in Section 2, it can be noticed that the calculated points in Figure 2, whose position will be advanced to the next time step, are located at the center of the shoreline segment, not at the two edges of the control volume.
Once the center nodes of the shoreline segment at the next time step are calculated, the control volume should be re-constructed.From the descriptions of the proposed evolution formulation and the numerical procedures, it is evident that the proposed model can easily deal with problems with discontinuities of the geometry.
Water 2022, 14, 3504 9 of 22 calculated points in Figure 2, whose position will be advanced to the next time step, located at the center of the shoreline segment, not at the two edges of the control volum Once the center nodes of the shoreline segment at the next time step are calculated, control volume should be re-constructed.From the descriptions of the proposed evolut formulation and the numerical procedures, it is evident that the proposed model can e ily deal with problems with discontinuities of the geometry.In this test, the wave height of normal incidence is 1m and we adopted the EEM [29,30], which can account for a rapidly varying topography and wave energy dissipat in the surf zone, to calculate plane wave fields.the parameters are adopted in the follo  In this test, the wave height of normal incidence is 1m and we adopted the EEMSE [29,30], which can account for a rapidly varying topography and wave energy dissipation in the surf zone, to calculate plane wave fields.the parameters are adopted in the following test: 100 , and λ = 0.43.In the first test, three different time increments, (∆t = 1, ∆t = 0.1, ∆t = 0.01, s), are used when m = 109 is adopted.The numerical results at t = 50 years are presented in Figure 6a and these three profiles of shoreline planform are almost identical to each other.The numerical comparison, provided in Figure 6a, can evidently verify the stability of the proposed model.Furthermore, the similar numerical experiments are implemented for m = 179 and m = 229, as well as the numerical comparisons are shown in Figure 6b,c.The excellent stability of the proposed numerical model is verified again in these two figures.In Figure 6a-c, it is evident that the profiles of steady shoreline are very different from the initial shorelines and the zones of advance and retreat can be simply identified.Moreover, the steady numerical results are compared well with each other in Figure 7

Numerical Validation and Comparisons
In this section, three cases are provided to show the numerical validation and engineering application of the proposed evolution formulation.In the first case, the numerical results are compared with previous study to demonstrate the effectiveness of the proposed model.In the second case, the total area of shoreline retreats and advances will be examined so as to prove the validation of conservation law of mass in the proposed model.Finally, the proposed evolution formulation is adopted to predict the change of shoreline profile in Yuguang Island Beach.The details of these three cases are clearly described in the following subsections.

Case 1
In order to validate the effectiveness of the proposed model, we used an example (Case B) in [20] as the Case 1 in this paper.The parameters of this experiment is presented in Table 1.In [20], they compared the numerical results with experimental data from [31].The shoreline profile at time T = 14 hr is adopted as the initial condition in this experiment, while the simulation is terminated at T = 38 hr.By using the proposed evolution formulation of Equation ( 7), the movement of shoreline planform from T = 14 hr to T = 38 hr can be accurately and efficiently calculated.The numerical results of shoreline profiles are depicted in Figure 8.In addition, the numerical solutions of [20] and the experiment data of [28] are also presented in the same figure.From the comparisons in Figure 8, it is evident that the numerical results by using the proposed evolution formulation are very similar with numerical solution in [20] and experimental data in [28].The effectiveness of

Numerical Validation and Comparisons
In this section, three cases are provided to show the numerical validation and engineering application of the proposed evolution formulation.In the first case, the numerical results are compared with previous study to demonstrate the effectiveness of the proposed model.In the second case, the total area of shoreline retreats and advances will be examined so as to prove the validation of conservation law of mass in the proposed model.Finally, the proposed evolution formulation is adopted to predict the change of shoreline profile in Yuguang Island Beach.The details of these three cases are clearly described in the following subsections.

Case 1
In order to validate the effectiveness of the proposed model, we used an example (Case B) in [20] as the Case 1 in this paper.The parameters of this experiment is presented in Table 1.In [20], they compared the numerical results with experimental data from [31].The shoreline profile at time T = 14 h is adopted as the initial condition in this experiment, while the simulation is terminated at T = 38 h.By using the proposed evolution formulation of Equation ( 7), the movement of shoreline planform from T = 14 h to T = 38 h can be accurately and efficiently calculated.The numerical results of shoreline profiles are depicted in Figure 8.In addition, the numerical solutions of [20] and the experiment data of [28] are also presented in the same figure.From the comparisons in Figure 8, it is evident that the numerical results by using the proposed evolution formulation are very similar with numerical solution in [20] and experimental data in [28].The effectiveness of the proposed model is verified from these comparisons in Figure 8.

Case 2
The Case 2 in Section 4.2 is the same as the numerical example in Section 3 in order to verify the accuracy of numerical results in Section 3. The accuracy of numerical results obtained by the proposed evolution formulation are carefully examined in this subsection, since the consistency and the stability of the proposed model are verified in Section 3. Two of the three straight lines in the initial shoreline of the numerical example, which is presented in Figure 5, are parallel to the y axis, so it is extremely difficult to apply the GEN-ESIS to simulate the dynamic transition of the sandy shoreline.Therefore, in the section, we used the hybrid formulation [20] to validate the accuracy of the proposed model.When the hybrid formulation is adopted to solve this example, the whole shoreline should be divided into two hooked zones and one unhooked zone, which are depicted in Figure 5, due to the appearance of these two breakwaters.Thus, two interfaces between different coordinates should be treated for the hybrid formulation.In the example, m = 229 and ∆t = 0.01 (sec) are adopted.The steady shorelines at t =50 years by using the proposed method and the hybrid formulation [20] are compared with each other in Figure 9.These two numerical solutions in Figure 9 are almost identical to each other, which can evidently validate the accuracy of the proposed model.It is interesting and reasonable to notice that the left half shoreline and the right half shoreline in Figure 9 are symmetric since all conditions in this example are symmetric.It is obviously verified that to place the origin of the polar coordinate at the tip of up-coast breakwater will not lead to unsymmetrical solutions.

Case 2
The Case 2 in Section 4.2 is the same as the numerical example in Section 3 in order to verify the accuracy of numerical results in Section 3. The accuracy of numerical results obtained by the proposed evolution formulation are carefully examined in this subsection, since the consistency and the stability of the proposed model are verified in Section 3. Two of the three straight lines in the initial shoreline of the numerical example, which is presented in Figure 5, are parallel to the y axis, so it is extremely difficult to apply the GENESIS to simulate the dynamic transition of the sandy shoreline.Therefore, in the section, we used the hybrid formulation [20] to validate the accuracy of the proposed model.When the hybrid formulation is adopted to solve this example, the whole shoreline should be divided into two hooked zones and one unhooked zone, which are depicted in Figure 5, due to the appearance of these two breakwaters.Thus, two interfaces between different coordinates should be treated for the hybrid formulation.In the example, m = 229 and ∆t = 0.01 (s) are adopted.The steady shorelines at t = 50 years by using the proposed method and the hybrid formulation [20] are compared with each other in Figure 9.These two numerical solutions in Figure 9 are almost identical to each other, which can evidently validate the accuracy of the proposed model.It is interesting and reasonable to notice that the left half shoreline and the right half shoreline in Figure 9 are symmetric since all conditions in this example are symmetric.It is obviously verified that to place the origin of the polar coordinate at the tip of up-coast breakwater will not lead to unsymmetrical solutions.
In addition to the steady solution of shoreline planform, the proposed evolution formulation can accurately acquire the evolutionary process of shorelines from the initial condition to the steady solutions.In this verification of the accuracy of solution during the evolutionary process of shorelines, m = 179 and ∆t = 0.01 s are used.In addition, the numerical results of the proposed model at different time in Figure 10 are in excellent agreements with the numerical solutions by using of the hybrid formulation [20].In Figure 10, it can be easily noticed that the profiles of shoreline are gradually changed from horizontal to half-moon-bay shape.Comparing the overall trend of the shoreline within hooked zone and unhooked zone, the unhooked zone withstands more breaking wave energy, which drives the longshore sediment transport toward to hooked zone and results in shoreline retreat in unhooked zone and shoreline advance in hooked zone.Furthermore, it could be noticed that in this Case the initial shoreline is constructed by three straight lines.Therefore, there will be two corners, which are the discontinuities of the geometry, in the initial shoreline.While the proposed novel evolution formulation for dynamic shoreline planform is derived in Section 2, it can be noticed that the calculated points in Figure 2, whose position will be advanced to the next time step, are located at the center of the shoreline segment, not at the two edges of the control volume.Once the center nodes of the shoreline segment at the next time step are calculated, the control volume will be re-constructed.From the numerical results in Figure 10, the profiles at the two sharp corners are smoothed by the proposed numerical scheme while the wave energy reformed the shape of shoreline.From the numerical results and comparisons in Figures 9 and 10, it can be emphasized that the proposed evolution formulation can accurately calculate the temporal variations of shoreline in addition to the steady shoreline planform.The profile s of steady shoreline s by using the proposed mode l and the hybrid formulation [20].
In addition to the steady solution of shoreline planform, the proposed evolution formulation can accurately acquire the evolutionary process of shorelines from the initial condition to the steady solutions.In this verification of the accuracy of solution during the evolutionary process of shorelines, m = 179 and ∆t = 0.01 sec are used.In addition, the numerical results of the proposed model at different time in Figure 10 are in excellent agreements with the numerical solutions by using of the hybrid formulation [20].In Figures 10, it can be easily noticed that the profiles of shoreline are gradually changed from horizontal to half-moon-bay shape.Comparing the overall trend of the shoreline within hooked zone and unhooked zone, the unhooked zone withstands more breaking wave energy, which drives the longshore sediment transport toward to hooked zone and results in shoreline retreat in unhooked zone and shoreline advance in hooked zone.Furthermore, it could be noticed that in this Case the initial shoreline is constructed by three straight lines.Therefore, there will be two corners, which are the discontinuities of the geometry, in the initial shoreline.While the proposed novel evolution formulation for dynamic shoreline planform is derived in Section 2, it can be noticed that the calculated points in Figure 2, whose position will be advanced to the next time step, are located at the center of the shoreline segment, not at the two edges of the control volume.Once the center nodes of the shoreline segment at the next time step are calculated, the control volume will be re-constructed.From the numerical results in Figure 10, the profiles at the two sharp corners are smoothed by the proposed numerical scheme while the wave energy reformed the shape of shoreline.From the numerical results and comparisons in Figures 9 and 10, it can be emphasized that the proposed evolution formulation can accurately calculate the temporal variations of shoreline in addition to the steady shoreline planform.The accuracy of the proposed model during the evolutionary process and at the steady state are validated in comparison with the solutions of the hybrid formulation in Figures 9 and 10.It is obvious that the shoreline profile of the initial shoreline is very different to the shoreline profile of the steady shoreline.Both of the retreat and advance can be discovered by comparing the shorelines at initial time with at steady state.Since The accuracy of the proposed model during the evolutionary process and at the steady state are validated in comparison with the solutions of the hybrid formulation in Figures 9 and 10.It is obvious that the shoreline profile of the initial shoreline is very different to the shoreline profile of the steady shoreline.Both of the retreat and advance can be discovered by comparing the shorelines at initial time with at steady state.Since only the longshore sediment transport is considered in the proposed formulation and no longshore sediment transport is imposed at the two ends of shoreline, it can be expected that the amount volume of erosion should be equal to the amount volume of deposition along the entire shoreline.In order to check the balance between retreat and advance, we calculated both zone, which are enclosed by the initial shoreline and the steady shoreline, as well as these numerical results are tabulated in Table 2. Three different time increments (∆t = 1, 0.1, 0.01, s) and three different numbers of control volume are used in this comparison.In each test, the area of retreat and advance are calculated and listed in Table 2.The difference between them are also calculated by the following equation, We used positive value to denote the area of advance and negative value to denote the area of retreat, so the smaller ε means that the numerical model can achieve better conservation law of mass for sediment.The numerical results and comparisons are tabulated in Table 2.The numerical results of the present model by using different numbers of control volume are very satisfactory and stable.Furthermore, the numerical results by using different time increments are very similar to each, which can verify the stability of the proposed method.In addition, when m = 229 is used, ε by using the proposed scheme is less than ε by using the hybrid scheme, which means that the proposed formulation outperforms the hybrid formulation in term of the conservation of sediment.
In the following test, we adopted three different angles of the incident wave, θ int , while m = 179 and ∆t = 0.01 s are used.The profiles of steady shorelines are demonstrated in Figure 11 for θ int = 0 • , ±30 • .The angle of incident wave implies the direction of wave energy, so the profiles of steady shorelines by using different angles of incident wave will be different from each other.In the numerical results of Figure 11, it is obvious that the angle of incident wave has great influence on the dynamic movements of sandy shoreline.The part of shoreline, which is directly attacked by incident wave without any shelter, will more toward offshore side.In the meantime, the other part of shoreline will be moved toward the land side.This phenomenon of retreats and advances by adopting different θ int is obvious in Figure 11 and has ever described by [32].Besides, it is interesting that the profiles of shoreline by using θ int = 30 • and θ int = −30 • are anti-symmetric.
strated in Figure 11 for   = 0 , ±30 .The angle of incident wave implies the direction of wave energy, so the profiles of steady shorelines by using different angles of incident wave will be different from each other.In the numerical results of Figure 11, it is obvious that the angle of incident wave has great influence on the dynamic movements of sandy shoreline.The part of shoreline, which is directly attacked by incident wave without any shelter, will more toward offshore side.In the meantime, the other part of shoreline will be moved toward the land side.This phenomenon of retreats and advances by adopting different   is obvious in Figure 11 and has ever described by [32].Besides, it is interesting that the profiles of shoreline by using   = 30 ° and   = −30 ° are anti-symmetric.In the next experiment, the angles of incident wave are gradually changed so as to realize the influence of the angles of incident wave to the shoreline planform in this test.  = −30 °,   = −20 °,   = −10 ° and    = 0 ° are adopted and the steady-state numerical results are shown in Figure 12.The profile of shoreline is symmetric with respect to the central vertical line for   = 0 °.When we gradually increase the angles of incident wave, it is apparent that the direction of wave energy shifted and the profiles of shoreline become anti-symmetric.As a result, the influence of angle of incident wave to the profile of shoreline can be accurately simulated by the proposed evolution formulation.In the next experiment, the angles of incident wave are gradually changed so as to realize the influence of the angles of incident wave to the shoreline planform in this test.
• and θ int = 0 • are adopted and the steady-state numerical results are shown in Figure 12.The profile of shoreline is symmetric with respect to the central vertical line for θ int = 0 • .When we gradually increase the angles of incident wave, it is apparent that the direction of wave energy shifted and the profiles of shoreline become anti-symmetric.As a result, the influence of angle of incident wave to the profile of shoreline can be accurately simulated by the proposed evolution formulation.
In Figure 13, we used a 1 ∼ a 5 to denote five zones, which are enclosed by the initial shoreline and the steady-state shoreline.From the numerical results in Figure 12, it is obvious that a 1 , a 3 , and a 5 are the zones of retreats, as well as a 2 , and a 4 are the zones of advance.In Table 3, the areas of a 1 ∼ a 5 are tabulated while different angles of incident wave are adopted.We noticed that the simulations need more time to achieve steady state if the angle of incident wave is larger than normal incidence.
In addition, it is natural to expect that a 1 = a 5 and a 2 = a 4 if the incident wave is normal to the shoreline.It is obvious that the imbalance of numerical results, a 1 = a 5 and a 2 = a 4 , appears in Table 3 for the case of normal incidence.There are two possible causes of these imbalance.One of these reasons is that the proposed evolution formulation may accumulate numerical errors for shoreline segments, which are far away from the origin of the coordinate system, since the r direction is not identical to the normal direction of boundary segment.The other possibility is the numerical error for calculating these five areas.Since in one-line model we only have few distinct spatial positions to represent the entire shoreline, it can be expected that the numerical quadrature for area may include apparent errors.Since the number of nodes along the entire shoreline is very less, the main reason for the imbalance of individual area (a1 = a5, a2 = a4) may be the imprecise numerical quadrature for area.Although the imbalance of individual area can be easily found in Table 3, the relative errors by Equation ( 11) is less than 3%, which means that the provided numerical results in this paper are acceptable and satisfactory.In Figure 13, we used  1 ~5 to denote five zones, which are enclosed by the initial shoreline and the steady-state shoreline.From the numerical results in Figure 12, it is obvious that  1 ,  3 , and  5 are the zones of retreats, as well as  2 , and  4 are the zones of advance.In Table 3, the areas of  1 ~5 are tabulated while different angles of incident wave are adopted.We noticed that the simulations need more time to achieve steady state if the angle of incident wave is larger than normal incidence.
In addition, it is natural to expect that  1 =  5 and  2 =  4 if the incident wave is normal to the shoreline.It is obvious that the imbalance of numerical results,  1 ≠  5 and  2 ≠  4 , appears in Table 3 for the case of normal incidence.There are two possible causes of these imbalance.One of these reasons is that the proposed evolution formulation may accumulate numerical errors for shoreline segments, which are far away from the origin of the coordinate system, since the r direction is not identical to the normal direction of boundary segment.The other possibility is the numerical error for calculating these five areas.Since in one-line model we only have few distinct spatial positions to represent the entire shoreline, it can be expected that the numerical quadrature for area may include apparent errors.Since the number of nodes along the entire shoreline is very less, the main reason for the imbalance of individual area (a1 ≠ a5, a2 ≠ a4) may be the imprecise numerical quadrature for area.Although the imbalance of individual area can be easily found in Table 3, the relative errors by Equation ( 11) is less than 3%, which means that the provided numerical results in this paper are acceptable and satisfactory.The time span, that is required for   = −30 ° to reach steady state, is adopted for all the numerical tests in Table 3.It is interesting to find that the area of  2 increased and the area of  4 decreased by increasing the angle of incident wave, though both of  2 and  4 are zones of advance.In addition, the absolute values of area of  1 and  5 increased and the absolute values of area of  3 decreased while the angle of incident wave is gradually enlarged.Furthermore, the conservation of sediment is also examined in Table 3.The excellent balance between retreats and advance, ε, are presented in this table.The effectiveness of the proposed evolution formulation is verified since ε is very small in comparing with the total amount of area of retreats or advance.The time span, that is required for θ int = −30 • to reach steady state, is adopted for all the numerical tests in Table 3.It is interesting to find that the area of a 2 increased and the area of a 4 decreased by increasing the angle of incident wave, though both of a 2 and a 4 are zones of advance.In addition, the absolute values of area of a 1 and a 5 increased and the absolute values of area of a 3 decreased while the angle of incident wave is gradually enlarged.Furthermore, the conservation of sediment is also examined in Table 3.The excellent balance between retreats and advance, ε, are presented in this table.
The effectiveness of the proposed evolution formulation is verified since ε is very small in comparing with the total amount of area of retreats or advance.
In the previous discussions, the five identical areas between the initial shoreline and the steady shorelines can be easily noticed by observing the numerical results in Figure 10.Though the satisfactory numerical results and comparisons are tabulated in Table 3, we adopted different lengths of initial shoreline, (L, S) to further examine the merits of the proposed evolutionary formulation.In Figure 14, two plane areas are defined.One is A 1 , which is the area enclosed by the initial shoreline and these two headlands.The other one is A 2 , which represents the area enclosed by the steady-state shoreline and the headlands.Although the shapes of A 1 and A 2 are distinct from each other, the values of A 1 and A 2 should be exactly the same, which means that the total areas of shoreline advances should be identical to the whole areas of shoreline retreats.In this test, three different m and three different ∆t are adopted to ensure the consistency and the stability, while three L and three S are used to form various aspect-ratio initial shorelines.The resultant numerical results are tabulated in Table 4.The relative error of total area, , is defined by the following equations, From the results, presented in Table 4, for the cases of fixed L and fixed S, the merical errors will decrease while m is increased, Meanwhile, the numerical errors decrease as ∆t is decreased for the cases of fixed L and fixed S. The excellent consiste and stability of the proposed formulation is evidently validated from these results.thermore, it is interesting to find that to increase either L or S will result in larger merical errors for the cases of fixed m and fixed ∆t.To increase either L or S wil crease the total length of initial shoreline and the representative length of each con volume, which is shown in Figure 3, when we used a fixed value of m for spatial disc zation.Thus, from the results in Table 4 it is reasonable to notice that the numerical e is proportional to the representative length of control volume.In addition, by compa the numerical result in Table 4 with Table 3, it is sensible to deduce that the large er of five areas and the imbalance of numerical results in Table 3 are mainly resulted f less nodes for numerical quadrature.From the results, presented in Table 4, for the cases of fixed L and fixed S, the numerical errors will decrease while m is increased, Meanwhile, the numerical errors will decrease as ∆t is decreased for the cases of fixed L and fixed S. The excellent consistency and stability of the proposed formulation is evidently validated from these results.Furthermore, it is interesting to find that to increase either L or S will result in larger numerical errors for the cases of fixed m and fixed ∆t.To increase either L or S will increase the total length of initial shoreline and the representative length of each control volume, which is shown in Figure 3, when we used a fixed value of m for spatial discretization.Thus, from the results in Table 4 it is reasonable to notice that the numerical error is proportional to the representative length of control volume.In addition, by comparing the numerical result in Table 4 with Table 3, it is sensible to deduce that the large errors of five areas and the imbalance of numerical results in Table 3 are mainly resulted from less nodes for numerical quadrature.

Case 3
In Case 3, the proposed evolution formulation is applied to an engineering application in Taiwan.In the south-west part of Taiwan, Yuguang Island Beach is located between Anping and Anping New Port of Tainan.It was found in google earth that the revetment in this section of coastline was not built in 1984.At the same time, Anping and Anping New Port were not built yet.In Figure 15, we could observe the straight shoreline of Yuguang Island Beach (red line) at 1985 was a parallel coastline.The two offshore breakwater (blue line) was built at 2000.Due to the existence of the offshore breakwater, after 20 years shoreline was advancing seaward (pink dot line).Based on the history of the above shoreline changes at Yuguang Island Beach, we know this area had raised beach nourishment in "Anping Port Adjacent Coastal Improvement Planning Report", by the picture it can be found that the offshore breakwater can effectively accumulate sediment and store sand in the sheltered areas behind the offshore breakwater.

Case 3
In Case 3, the proposed evolution formulation is applied to an engineering application in Taiwan.In the south-west part of Taiwan, Yuguang Island Beach is located between Anping and Anping New Port of Tainan.It was found in google earth that the revetment in this section of coastline was not built in 1984.At the same time, Anping and Anping New Port were not built yet.In Figure 15, we could observe the straight shoreline of Yuguang Island Beach (red line) at 1985 was a parallel coastline.The two offshore breakwater (blue line) was built at 2000.Due to the existence of the offshore breakwater, after 20 years shoreline was advancing seaward (pink dot line).Based on the history of the above shoreline changes at Yuguang Island Beach, we know this area had raised beach nourishment in "Anping Port Adjacent Coastal Improvement Planning Report", by the picture it can be found that the offshore breakwater can effectively accumulate sediment and store sand in the sheltered areas behind the offshore breakwater.Before the numerical simulation, the wave data in the considered region over the past years were collected from the Central Weather Bureau of Taiwan and these data are presented in Table 5.By averaging these wave data over 20 years, it is known that the average wave height is equal to 0.66 m, wave period is 5.5 s, and wave direction is SW.Although the hourly wave data in the study site are available, we still used the average wave height due to the following reason.In Case 3, some engineering projects of beach nourishment have been implemented in the study site during the past decades, but the details of the beach nourishment are unavailable.Thus, we used the cross-shore sediment transport, q, to replace the effectiveness of beach nourishment.According to the basic assumption of the proposed formula, the volume in the calculation domain is conserved, since no source of sand is found from offshore area and two ends of shoreline.Obviously, the volume of the beach in the case is not conserved and there should be some source of sand from the offshore area.Therefore, the supplementary term q of the external source of clothing needs to be added to the proposed formulation, the supplementary term q in the Case 3 is calculated as 0.00166 m 2 /h, and the proposed evolution formulation is rewritten as follows, In the numerical simulation, the shoreline in 1985 was used as the initial shoreline, and the simulation time was 20 years.The simulation results are shown by the yellow line in Figure 15.From the result in Figure 15, it can be found that the shoreline changed from a parallel coastline to a half-moon shoreline, which is similar to the result in Section 4.2.In Figure 15, the numerical result of shoreline profile is in good agreement with real shoreline position.The numerical results and comparisons in Case 3 reveals that the proposed evolution formulation can be efficiently applied to real engineering applications.It is unreal to use a constant supplementary term q in real application, so the proper description of supplementary term q will be studied in our future research.From the numerical results and comparisons in Case 3, it can be found that to use cross-shore sediment transport to replace the amount of beach nourishment is effective.In addition, the proposed novel Water 2022, 14, 3504 20 of 22 formulation can include both of the cross-shore sediment transport and longshore sediment transport at the same time.

Conclusions
In this paper, a novel evolution formulation, based on the polar coordinate, is proposed to accurately simulate the dynamic shoreline planform of crenulate-shaped bay.The integrated nature forces of wind, wave, current, tide et al. will continuously result in shoreline retreats or advances along sandy beach, so the proposed evolution formulation to forecast the dynamic movement of shoreline plays an important role in engineering design and coastal protection.The proposed evolution formulation is derived by considering the longshore sediment transport and the rectangular control volume.A simple timemarching scheme is presented in the proposed formulation due to the use of the rectangular control volume.Furthermore, the entire shoreline can be described by using only the polar coordinate, which can simplify the numerical implementations in comparison with other formulations.Consequently, a simple, accurate and efficient numerical scheme is proposed in this paper to simulate the dynamic movements of sandy shoreline.
In order to verify the merits of the proposed numerical scheme, different time increments and different numbers of control volume are adopted in numerical simulation.The numerical comparisons of results between the proposed scheme with other formulations are satisfactory and can be used to validate the accuracy of the proposed method.Moreover, the total amount of area of retreat and advance are calculated in order to demonstrate the conservation of sediments in the proposed formulation.Besides, different offshore wave incidence directions are examined by the proposed method, and the influence of the angles of incident wave to the shoreline planform is clearly revealed by numerical experiments.In addition to the comparison between the numerical and the laboratory data, we used a field case in Yuguang Island Beach to show the feasibility of the proposed model to engineering application.
From the numerical comparisons, it is verified that the proposed model can well represent the changes of the shoreline under the long-term action of waves.However, there are still some placed of the proposed evolution formulation can be improved, such as the control volume deformation caused by the distance from the calculation point in the downcoast, and uneven distribution of calculated points on the shoreline.These improvements and modifications will be the primary focus of our follow-up research.Since the mathematical derivation of the stability condition of the proposed formulation is very essential to engineering applications, it is an ongoing subject of our research the mathematical derivation of stability criterion of the proposed formulation.

Figure 1 .
Figure 1.The sche matic diagram of the shoreline planform of the cre nulate-shaped bay.

Figure 1 .
Figure 1.The schematic diagram of the shoreline planform of the crenulate-shaped bay.Wa ter 2022, 14, x FOR PEER REVIEW 6 of 22

Figure 2 .
Figure 2. The control volume be twe en two consecutive time steps.

Figure 2 .
Figure 2. The control volume between two consecutive time steps.

Figure 3 .
Figure 3.The re ctangular control volume de fine d in this study and its dime nsion.

Figure 3 .
Figure 3.The rectangular control volume defined in this study and its dimension.

Figure 4 .
Figure 4.The sche matic diagram for diffe rent angles used in this study.

Figure 4 .
Figure 4.The schematic diagram for different angles used in this study.

1 100,
ing test:  1 = 0.77,  2 = 0.5 1 ,   = 10,  = 1033 ,   = 2650 ,  =   = 0.43 the first test, three different time increments, (Δt = 1, Δt = 0.1, Δt = 0.01, sec), are us when  = 109 is adopted.The numerical results at t = 50 years are presented in Figu 6a and these three profiles of shoreline planform are almost identical to each other.T numerical comparison, provided in Figure 6a, can evidently verify the stability of the p posed model.Furthermore, the similar numerical experiments are implemented for  179 and  = 229 , as well as the numerical comparisons are shown in Figure 6b,c.T excellent stability of the proposed numerical model is verified again in these two figur In Figure 6a-c, it is evident that the profiles of steady shoreline are very different from initial shorelines and the zones of advance and retreat can be simply identified.Moreov
when m = 109, m = 179 and m = 229 are adopted.The superior consistency of the proposed model is clearly verified from these figures, since the solutions by using different numbers of control volume are in good agreements with each other.Consequently, the merits of the proposed model, including the consistency and the stability, are clearly verified by numerical experiments.thesteady numerical results are compared well with each other in Figure7when  = 109,  = 179 and  = 229 are adopted.The superior consistency of the proposed model is clearly verified from these figures, since the solutions by using different numbers of control volume are in good agreements with each other.Consequently, the merits of the proposed model, including the consistency and the stability, are clearly verified by numerical experiments.

Figure 6 .
Figure 6.The numerical results and comparison by using different ∆t for (a) m = 109, (b) m = 179 and (c) m = 229.

Figure 6 .
Figure 6.The numerical results and comparison by using different ∆t for (a) m = 109, (b) m = 179 and (c) m = 229.

Figure 7 .
Figure 7.The profile s of initial shore line and steady shoreline by using diffe rent numbers of control volume (∆t = 0.01 sec).

Figure 7 .
Figure 7.The profiles of initial shoreline and steady shoreline by using different numbers of control volume (∆t = 0.01 s).

22 Figure 9 .
Figure 9.The profile s of steady shoreline s by using the proposed mode l and the hybrid formulation [20].

Figure 9 . 22 Figure 10 .
Figure 9.The profiles of steady shorelines by using the proposed model and the hybrid formulation [20].Wa ter 2022, 14, x FOR PEER REVIEW 14 of 22

Figure 11 .
Figure 11.The profile of ste ady-state shoreline s by adopting three diffe rent angles of incide nt wave.(m = 179, ∆t = 0.01 sec).

Figure 11 .
Figure 11.The profile of steady-state shorelines by adopting three different angles of incident wave.(m = 179, ∆t = 0.01 s).

Figure 12 .
Figure 12.The profile s of ste ady-state shoreline by using four diffe re nt angles of incide nt wave.(m = 179, ∆t = 0.01 sec).

Figure 13 .
Figure 13.The schematic diagram for the five zones of shoreline retreats and advance.

)Figure 14 .
Figure 14.The de finitions of two distinct are as,  1 and  2 , which are e nclosed by shorelines two he adlands.

Figure
Figure The definitions of two distinct areas, A 1 and A 2 , which are enclosed by shorelines and two headlands.

Figure 15 .
Figure 15.The profile s of shore line in Yuguang Island Be ach at 1985 and 2021.Figure 15.The profiles of shoreline in Yuguang Island Beach at 1985 and 2021.

Figure 15 .
Figure 15.The profile s of shore line in Yuguang Island Be ach at 1985 and 2021.Figure 15.The profiles of shoreline in Yuguang Island Beach at 1985 and 2021.

Table 1 .
The physical parameters of Case 1.

Table 1 .
The physical parameters of Case 1.

Table 2 .
Numerical results of area of shoreline retreats and advance by comparing the initial shoreline with steady shoreline.

Table 3 .
The amount area of shoreline retreats and advance for different angles of incident waves.

Table 3 .
The amount area of shoreline retreats and advance for different angles of incident waves.

Table 4 .
The numerical error of areas by using different m, ∆t, S and L.

Table 5 .
Parameters and results of Case 3.