1 INTERACTION BETWEEN TWO CIRCULAR SANDY ISLANDS ON FLAT SHALLOW SEABED OWING TO WAVES

On a flat shallow seabed, sand spits and cuspate forelands with rhythmic shapes may develop and a barrier island can elongate alongshore, which significantly differs from the beach changes on coasts facing a deep ocean. In this study, the interaction between two circular sandy islands on a flat shallow seabed owing to waves was investigated by numerical simulation, focusing on the wave-sheltering effect of the islands themselves. The topographic changes caused by the interaction between two sandy islands were predicted using the BG model (a three-dimensional model for predicting beach changes based on Bagnold’s concept).


INTRODUCTION
showed that the formation of Nantasket Beach in Massachusetts can be explained by the process whereby an island originally formed by the deposition of drumlin in a shallow sea was eroded by waves, and bed materials supplied from this erosion were transported alongshore forming a slender sand bar, resulting in the connection of the island to another island (Fig. 1).This example shows that a sand bar may extend alongshore in a shallow water body, and also two small islands may connect with each other by the extension of a comet tail behind the island.Although Miyahara et al. (2013) investigated the topographic changes around a circular island surrounded by a deep sea, characteristic beach changes, which significantly differ from those observed under such a condition, can often be observed in a shallow water body: the development of sand spits and cuspate forelands with rhythmic shapes, as shown by Serizawa et al. (2012) and Uda et al. (2014), and the elongation of a barrier island (Serizawa and Uda, 2011).However, the deformation and interaction of islands composed of sand in a shallow sea have not yet been studied, and the understanding of the wave-sheltering effect by the islands themselves is still insufficient.In this study, we investigate the deformation and interaction of multiple islands composed of sand in a shallow sea, as the example given by Davis and Fitzgerald (2004) using the BG model proposed by San-nami et al. (2013), which is a three-dimensional (3-D) model for predicting beach changes based on Bagnold's concept.

CONNECTED ISLANDS IN HINGHAM BAY IN MASSACHUSETTS
Figure 2 shows the satellite image of Hingham Bay near Boston Harbor, where Davis and Fitzgerald (2004)

MODEL FOR PREDICTING BEACH CHANGES
The sand transport equation employed in this study is Eq.(1), using the expression of the wave energy evaluated at the breaking point, the same as in the BG model proposed by Serizawa et al. (2006) and San-nami et al. (2013).The variables in Eq. ( 1) are given by Eqs. ( 2) -( 10). tan I ε i+1 ( ) = min Eq.( 6) , Here, q = q x , q y ( ) is the net sand transport flux, Z (x, y, t) is the seabed elevation with reference to the still water level (Z = 0), ∇Z = ∂Z ∂ x, ∂Z ∂ y ( ) is the seabed slope vector, e w is the unit vector of wave direction, α is the angle between the wave direction and the direction normal to the contour line, x w is the coordinate along the direction of wave propagation, tanβ w is the seabed slope measured along the direction of wave propagation, tanβ c is the equilibrium slope of sand, and K s is the longshore and crossshore sand transport coefficient.C 0 is the coefficient transforming the immersed weight expression to the volumetric expression ( } ; ρ is the seawater density, ρ s is the specific gravity of sand, p is the porosity of sand, and g is the acceleration due to gravity), h c is the depth of closure, and h R is the berm height.ε(Z) in Eq. ( 5) is the depth distribution of sand transport, (EC g ) b is the wave energy flux at the breaking point, H b is the breaker height, H 1/3 is the significant wave height of the incident waves, and γ is the ratio of the breaker height relative to the water depth.
The index i in Eqs. ( 8) and ( 9) is the mesh number along the x w -axis, and min means the selection of the smaller of either value in parentheses.In addition, k 1 = (4.004) 2 in Eq. ( 10b) is a constant in the relationship between the wave energy E and the significant wave height when the probability of the wave height of irregular waves is assumed to be given by the Rayleigh distribution (Horikawa, 1988).In the calculation, the local beach slope measured along the wave ray was used as the beach slope in Eq.
For the calculation of the P value, another coordinate system different from that for the calculation of beach changes was used (San-nami et al., 2013), in which the x w -and y w -axes were taken along the wave direction and the direction normal to that, respectively (Fig. 5).A fixed coordinate system of (x, y) is for the calculation of beach changes with a calculation domain of a rectangular ABCD, whereas the P value was calculated in the rectangular domain of A'B'C'D' including the domain ABCD with the coordinate system of (x w , y w ).In the wave calculation, wave refraction was neglected and waves were assumed to propagate in a straight manner while maintaining wave incident angle.
The distance along the x w -axis is subdivided by the mesh of Δx w , and a cumulative function of ε(Z) I ε (Z) is introduced as Eq. ( 6).Here, ε(Z) is assumed to have a uniform distribution (Eq.( 5)) and I ε (Z) is a function that takes unity in a zone deeper than the depth of closure, decreases with the water depth, and becomes 0 in a zone higher than the berm height.Using Eq. ( 6), Eq. ( 2) is transformed into Eq.( 7).Thus, the P value can be determined from the derivative of I ε (Z) at a point along the x w -axis using Eq. ( 7).I ε (Z) is calculated along the x w -axis from the starting point of wave incidence toward the direction of wave propagation; I ε (Z) at the (i+1) th point can be calculated from Eq. ( 8) with the given value of I ε (Z) at the i th point, when the initial value of I ε (Z) at the offshore end is given and the mesh location is denoted by x w i ( ) = iΔx w .
Furthermore, the P value at the (i+1/2) th point is calculated from Eq. ( 9) which is the derivative expression of Eq. ( 7).Note that in calculating I ε (Z) at the (i+1) th point, the smaller value between the one calculated from Eq. ( 6) given the elevation Z(i+1) at the (i+1) th point and I ε (Z) at the i th point is adopted.As a result, the value of I ε (Z) corresponding to the minimum water depth (maximum Z) between the offshore end to a designated point along the x w -axis can be adopted.
Using this procedure, the depth distribution along the x w -axis between the offshore end and a certain point is automatically taken into account in the calculation of the P value.For example, when there is a location with an elevation higher than the berm height, the P value in the shoreward zone is automatically reset as 0 with no regard to the elevation and water depth.Furthermore, this procedure becomes important when the beach profile along the x w -axis has several uneven shapes.Consider the case wherein the x w -axis passes through near the tip of the two sand spits.The wave energy is reduced when waves pass through near the tip of the first sand spit, resulting in the reduction of the wave energy to the second sand spit.Using this procedure, the wave-sheltering effect due to the existence of the multiple sand spits can be automatically evaluated in contrast to the method whereby the P value is calculated by substituting the local elevation Z into ε(Z) in Eq. (2).Furthermore, the P value on the impermeable breakwater is set 0 to take the wave-sheltering effect of the breakwater into account.
In addition, the P value integrated from a location on land where the elevation exceeds the berm height to an offshore end along the x w -axis is always equivalent to the wave energy flux at the breaking point (EC g ) b with no regard to the seabed topography.Because (EC g ) b corresponds to the entire power of the incident waves, the exact satisfaction of this condition is reasonable.
The calculation of the P value is independently carried out along each x w -axis.The P value calculated at a point of (x w , y w ) is memorized and the P value at a point of (x, y) necessary for the calculation of beach changes was interpolated from the memorized value at a point of (x w , y w ).The mesh intervals of Δx w and Δy w are the same as those of Δx and Δy.Here, Δx and Δy are the mesh intervals of the coordinate system for the calculation of beach changes, and we assumed the equivalent condition of Δx = Δy.The beach changes were calculated by explicitly solving the continuity equations ( ∂Z ∂t + ∇ • q = 0 ) on the staggered meshes using the sand transport fluxes obtained from Eq. ( 1).The P value, which is subject to change with the propagation of waves, was recalculated at every calculation step of beach changes.In addition, when the ground elevation Z approaches the upper limit (h R ) or the lower limit of beach changes (-h c ) in the calculation process of beach changes, the sand transport rate was reduced by multiplying the reduction ratio using the method of Uda et al. (2013), so that the beach changes over these limits do not occur.

CALCULATION CONDITIONS
The circular islands composed of sand were placed on a plane solid seabed with a water depth of 2 m and the subsequent changes under the wave action were predicted.Taking into consideration that the beach changes occur in extremely shallow water, the depth distribution of the sand transport was assumed to be uniform.Two wave incident conditions were considered: a unidirectional wave incidence and the wave incidence with the direction ranging between 0 and 360 degrees with the same probability.The incident wave height H 1/3 is assumed to be 1 m in both cases.The other calculation conditions are shown in Table 1.Boundary conditions Shoreward and landward ends q x = 0 Right and left boundaries q y = 0

RESULTS
Deformation of a Circular Island under Action of Unidirectional Waves Miyahara et al. (2013) predicted the topographic changes of a circular island composed of sand using the BG model.In their study, the deformation of an island with the formation of sand spits occurred in the seabed, the depth of which was sufficiently large compared with the depth of closure, h c , and therefore, sand was permitted to discharge into the deep sea.This resulted in the gradual retardation of the development of sand spits.In this study, to investigate the extension of a slender sand bar or sand spit from an island composed of sand in the shallow sea, as shown by Davis and Fitzgerald (2004), the deformation of a circular island under the action of the unidirectional waves was simulated first.
Consider a flat shallow seabed with a water depth of 2 m, which is equivalent to the depth of closure h c .An island composed of sand is assumed to be produced with an initial beach slope of 1/10 and the beach elevation of 1 m which is equivalent to the berm height h R .This flat shallow seabed is assumed to be a solid bed. Figure 6 shows the initial topography and the results of the calculation.
Owing to the wave action to the initial topography, the exposed side of the island was eroded, and a mildly curved shoreline was formed after 200 steps, and sand eroded from the exposed side was transported to both sides of the island, resulting in the formation of two couples of symmetric protrusions, A and B, as shown in Fig. 6(b).After 400 steps, protrusion B markedly extended, whereas protrusion A stopped extending because of the wave-sheltering effect due to sand spit B. The sand spits on both sides of the circular island continued to develop after 600 and 800 steps.After 1000 steps, the sand spits on both sides of the circular island further extended so that the wave-sheltering effect by those sand spits covered all the lees of the island.Thus, the exposed side of the circular island was eroded in an extremely shallow sea and sand produced from the erosion was laterally transported, resulting in the formation of oblique slender sand spits on both sides.This explains well the beach changes around the sandy island located in a shallow sea in Fig. 1, as shown by Davis and Fitzgerald (2004).

Deformation of a Circular Island behind an Impermeable Offshore Breakwater
Because the wave-sheltering effect is assumed to play a primary role in the deformation of the sandy island, a case in which an impermeable offshore breakwater was extended straightly near the circular island was first considered.Consider the condition that waves are incident to a circular sandy island from all the directions between 0 and 360 degrees with the same probability.Under these conditions, the topographic changes of a circular island were calculated, in which the wave direction was randomly determined every time step of the calculation of beach changes.
Consider a circular sandy island with the origin at a point of (x, y) = (0, 0) and the radius of 100 m, as shown in Fig. 7(a), and a straight breakwater of 300 m length is placed at x = -200 m.The elevation and slope of the island are assumed to be 1 m and 1/10, respectively, and the water depth surrounding the island is assumed to be 2 m.Although waves are incident from all directions between 0 and 360 degrees, the waves incident from the direction of the negative x-axis are sheltered by the breakwater, as schematically shown in Fig. 8, and the wave action from this direction to the sandy island is markedly reduced.On the other hand, waves incident from the direction of the positive x-axis are not affected by the breakwater, the same as the waves incident from the ± y-axis.Thus, gradual sand movement towards the breakwater is expected to occur along the shoreline of the circular island.When waves are incident to the initial topography shown in Fig. 7(a) under the above-mentioned conditions, a pair of small sand bars extended toward the offshore breakwater from the upper end of the island after 10 3 steps.After 2×10 3 steps, the sand bars formed between the island and the offshore breakwater further extended and the tip of the sand bars connected the offshore breakwater.After 5×10 3 steps, the opposite side of the island far from the breakwater was eroded, sand was transported toward the breakwater, and a pair of sand bars had fully connected to the breakwater leaving a small pond in the center of the sandy beach.With time elapsed as 1×10 4 and 2×10 4 steps, a large amount of sand was transported to the lee of the breakwater losing the shape of a circular island and, finally, after 5×10 4 steps, a sandy beach with of triangular shape was formed in the lee of the breakwater.

Interaction of Sandy Islands of Symmetric Form
Although in the case mentioned above the wave-sheltering effect was induced by the installation of an impermeable offshore breakwater, the same effect may be induced by the sandy islands themselves.We, therefore, consider two sandy islands with the origin at points of (x, y) = (0, 200) and (0, -200) and with the radius of 100 m, as shown in Fig. 9.It is assumed that the two islands have the same size, the elevation and beach slope of the islands are 1 m and 1/10, respectively, and the two islands exist in a shallow sea of 2 m depth, the same as in Case 2. Although the distance between the centers of the island is 400 m, the shorelines of the islands are 200 m apart.We assume that waves are incident to these islands from all the directions between 0 and 360 degrees with the same probability as in Case 2. Under these conditions, island A is subject to the wave-sheltering effect by island B when waves are incident from the direction of the positive y-axis, as schematically shown in Fig. 10.Similarly, island B is subject to the wave-sheltering effect by island A when waves are incident from the direction of the negative y-axis.On the other hand, when waves are incident from ±x-axis direction, the wave field around islands A and B is not subject to the wave-sheltering effect of either island.
Because waves are incident from all the directions around the islands, if waves are incident from the direction with a ±y component rather than the wave incidence from the direction normal to the yaxis, island A or B is subject to the wave-sheltering effect by island B or A, respectively.Since the offshore breakwater is located at a fixed position, the offshore breakwater creates a stationary waveshelter zone behind the structure, whereas in the case of multiple islands, the islands producing the wave-sheltering effects themselves may deform, resulting in a temporal change in the wave-sheltering effect.These are the major differences in both cases.
When waves were incident to the initial topography, as shown in Fig. 9(a), under the conditions mentioned above, slender sand bars started to extend from the islands in the direction opposite to each other between the islands after 10 3 steps.This development of a cuspate foreland in the lee of the island explains well the development of a sandy beach on the Maebama coast located in the western part of Miyako Island (Kikuchi et al., 2002) and the formation of Nagayama Beach on Irabu Island located west of Miyako Island (Uda et al., 2013).
After 2×10 3 steps, the sand bars extended from islands A and B almost connected with each other.After 5×10 3 steps, the two islands connected with each other and became a single island with a neck in the center.As time further elapsed at 1×10 4 and 2×10 4 steps, the neck in the center gradually disappeared and the island was reduced to an elliptic form.After 1×10 5 and 2×10 5 steps, an island of almost elliptic form was formed with a center at a point of (x, y) = (0, 0).

Interaction of Two Sandy Islands of Asymmetric Form
Although the islands of the same size produce the equivalent wave-sheltering effect, it significantly changes in the case of the islands of different sizes.Because the wave-sheltering effect of a large island is more effective than that of a small island, asymmetric beach changes may occur.Here, two circular sandy islands A and B both centered at points of (x, y) = (0, 200) and (0, -200), and with radii of 100 m (A) and 50 m (B) were considered, as shown in Fig. 11.In this case, the elevation of the islands and water depth surrounding the islands were assumed to be 1 m and 2 m, respectively.In addition, the beach slope is assumed to be 1/10, the same as in Case 3. The length between the centers of the islands is 400 m and the channel width between the islands is 250 m.Waves are incident from all the directions between 0 and 360 degrees to these islands with the same probability as in Case 3.Under these conditions, island A is subject to the wave-sheltering effect of island B when waves are incident from the direction of the positive y-axis, but the wave-sheltering effect is relatively weak because of the smaller scale of island B, as schematically shown in Fig. 12.On the other hand, under the wave incidence from the direction of the negative y-axis, island B is subject to the stronger wavesheltering effect of island A because island A is larger than island B. Since waves are incident from all the directions around the islands, if waves are incident from the direction with a ±y component, island A or B is subject to the wave-sheltering effect by island B or A, respectively, but the wave-sheltering effect always becomes asymmetric.
When waves are incident to the initial bathymetry (Fig. 11(a)) under the conditions mentioned above, the cuspate foreland started to extend in opposite directions between the islands after 10 3 steps, resulting in the larger cuspate foreland of island B to island A because island B is subject to the stronger wave-sheltering effect of island A. After 2×10 3 steps the sand bar markedly extended from island B to island A, and after 5×10 3 steps, the two islands extended and almost connected with each other.After 1×10 4 and 2×10 4 steps, a slender sand bar, the width of which decreases rightward, was formed.After a large number of time steps, such as 1×10 5 and 2×10 5 steps, the island became elliptical in shape, and the center of the island moved as a whole near the point of (x, y) = (-100, 0) leftward compared with that in the case that the scale of the islands is equivalent, as shown in Fig. 11(h).Thus, the wave-sheltering effect of the islands themselves plays an important role in the deformation of the islands located in shallow sea.Davis and Fitzgerald (2004) showed that the formation of Nantasket Beach in Massachusetts can be explained by the process whereby an island originally formed by the deposition of drumlin in a shallow sea was eroded by waves, and bed material supplied from this erosion was transported alongshore forming a slender sand bar, resulting in the connection of one island to another.This example shows that a sand bar may extend alongshore in a shallow water body, and also two small islands may connect with each other by the extension of a comet tail behind the island.In this study, the deformation and interaction of multiple islands composed of sand in a flat shallow sea were investigated using the BG model.The 3-D beach changes and interactions of these multiple islands could be predicted using this model, and the applicability of this model was further expanded.

Figure 1 .
Figure 1.Connecting process of islands formed by the deposition of drumlin in shallow sea (Davis and Fitzgerald, 2004).
Figure2shows the satellite image of Hingham Bay near Boston Harbor, whereDavis and Fitzgerald (2004) investigated the formation of Nantasket Beach.Selecting two rectangular areas 1 and

Figure. 2 .
Figure. 2. Satellite image of Hingham Bay in Massachusetts.

Figure
Figure. 5. Definition of coordinate system.

Figure. 6 .
Figure. 6. Deformation of a circular island composed of sand under the action of unidirectional waves.

Figure. 7 .
Figure. 7. Deformation of a circular island composed of sand behind an impermeable breakwater.

Figure
Figure.11.Deformation and interaction of two asymmetric islands composed of sand.