1. Introduction
Crenulate-shaped or headland bays are quite common on exposed sedimentary coasts containing headlands, and represent about 50% of the world’s coastline [
1,
2]. The term headland-bay beach has been used to define a shoreline bounded by rocky outcrops or headlands, either natural or man-made, which lead to diffraction of incoming waves. Particularly, the predominant waves are diffracted in such a way as to break simultaneously around the periphery of the bay once an equilibrium plan-shape has been established (static equilibrium condition, [
3]). The fame of the headland bay beaches, in fact, lies in their equilibrium condition, static or dynamic, which ensures that they are considered a way to achieve coastline stabilization [
4]. Static equilibrium is a condition characterized by the absence of littoral drift, without the need for sediment supply to preserve its long-term stability; on the other hand, a dynamic equilibrium condition requires sediment supply, from updrift and/or from another kind of source, to maintain its stability and not retreat towards the limit defined by static equilibrium position [
3].
Typically, the plan-shape of a single-headland bay is characterized by an upcoast curved zone (diffraction zone), a gentle transition zone and a relatively straight tangential segment on the downdrift end of the bay (illuminated zone), which is largely orthogonal to the dominant wave direction; this equilibrium plan-form is that assumed by the bay at a relatively long-term scale (e.g., annual to decadal) as a response to the predominant wave direction. Short-term fluctuations arising from beach-storm interactions, which could cause severe beach berm retreat, can be neglected due to their reversibility.
Devising a headland bay beach system in static equilibrium in such a way that the reorientation of the shoreline cancels the longshore sediment transport is called “
headland control”, and can be considered as a valuable option for coastal stabilization, which is frequently addressed via traditional detached low-crested breakwaters [
5,
6] or artificial reefs [
7,
8].
The headland control concept has been suggested for engineering use by [
9,
10] as a naturally functioning and preferable means of shore protection. In this regard, the formation of a crenulate-shaped bay on a sedimentary coastline, under oblique attack of persistent swell, is the most stable beach generated by nature [
11]. The headland control approach appears to be of the utmost importance in the field of coastal stabilization against beach erosion, which increasingly affects many parts of the world’s coastline. Hence, the employment of a crenulate-shaped bay to stabilization of a shoreline may be a powerful tool for engineering purposes.
Several studies and researches have been carried out in order to develop functional models describing static equilibrium bays’ shapes: the logarithmic-spiral model [
12], hyperbolic-tangent model [
13,
14] and parabolic model [
15,
16,
17]. However, they are merely of an empirical nature, lacking further insight on relationships between incident wave characteristics and beach shape. As a matter of fact, none of these shapes are derived directly from the acting physical processes that developed the shoreline; rather, they are observational. Consequently, despite the strength in the assessment of stability of existing beaches, they are affected by uncertainties in the design of new artificial beaches through the headland control approach, which instead requires a deep knowledge of the dominant physical processes that govern the plan-shape of the bay. Without that, the project design of a new beach could result extremely challenging.
With the aim of overcoming such drawbacks and, broadly, to establish a relationship between wave forcings (diffraction and refraction) and bay shape response, the main task of this paper is founding a possible correspondence between static equilibrium profiles and wave characteristics, accounting for the relationship between equilibrium shape profiles and wave fronts. In fact, it is commonly believed that equilibrium beach profiles follow the wave front trend, but this has not been proved in literature so far, and no research has clarified in depth how wave characteristics, and particularly wave fronts, could shape a crenulated stretch of coast. For these reasons, this research represents a first step towards a development of a guidance which could help in engineering and morphological practice.
The analysis has been carried out via numerical modeling, which has long been proved to be a powerful tool, suited to even complex hydrodynamic phenomena [
18,
19,
20].
Numerical experiments have been performed using the MIKE 21 Boussinesq Wave Module (BW) [
21], where wave fronts have been compared to the hyperbolic-tangent equilibrium profile, analysing the influence of wave direction, wave period and refraction phenomenon. A correspondence function, called the “wave-front-bay-shape equation” has been established, offering an easy application to engineering uses due to the simple geometric interpretation of its controlling parameters. It is noteworthy that, being based on wave front analysis, the employed approach is essentially linear (e.g., [
22]), although Boussinesq wave models account for weak non-linearity. As such, the parametric study discussed in
Section 3 treats the wave height as an invariant of the problem.
The simulation’s outcomes seem to indicate that equilibrium beach profiles of a single headland bay correspond to a simple translation of wave fronts normal to the propagation line at the headland tip. Moreover, the application of the “wave-front-bay-shape equation” to the case-study bay of the Bagnoli coast (south-west of Italy) is described in the paper.
2. Background
Since the beginning of nineteenth century, many coastal geologists, geographers and engineers have been trying to predict the shoreline plan-shape of headland-bay beaches. Despite the complexity involved in coastal processes, simple empirical expressions have been derived since the 1940s in order to fit part or whole of the bay periphery. Among these approaches, the logarithmic spiral model [
12], hyperbolic tangent model [
13,
14] and parabolic model [
15,
16,
17] have turned out to be the most acceptable expressions for practical applications to headland bay beaches in static equilibrium. However, none of these models are derived directly from the acting physical processes that developed the shape; rather, they are based on the observation of the shoreline plan-shape, and so are lacking in a correlation between shoreline response and wave forcings (refraction and diffraction).
The logarithmic-spiral model was introduced by [
12] observing an early 1940 imagery of Half-Moon Bay in California, USA. With field observations, the author realized that the bay adopted an equilibrium shape that is similar to a log-spiral. A definition sketch is given in
Figure 1. The log-spiral equation, in polar coordinates, is given by:
where:
Despite of being the very first model proposed, representing a milestone in the field of crenulated bays research, the log-spiral model has long been criticized [
3,
23] because of several drawbacks which make it tricky to use. First, due to its constant curvature, the equation does not fit the relatively straight section of the headland bay beach; secondly, the pole of the log-spiral curve must be found by trial and error, because it does not coincide with the diffraction point, and could even deviate from it by a distance ranging from centimetres to kilometres (see [
24]); therefore, it is not possible to predict the effect of relocating the headland (e.g., designing a headland bay beach introducing a coastal structure); and finally it does not take into account wave direction. In light of these disadvantages, two other empirical equations were developed: the hyperbolic-tangent equation and the parabolic equation.
The parabolic bay shape equation is a second-order polynomial equation developed by [
15] and [
16,
17] in two separate works, from fitting the plan-shape of 27 mixed cases of prototype and model bays believed to be in static equilibrium:
where:
= Control line length
= Radius to a point along the curve at an angle θ
= Wave obliquity
= Constants generated by regression analysis
= Angle between wave crest and radius to any point on the bay periphery
The two basic parameters are the reference wave obliquity
β and control line length
Rβ (
Figure 2). The variable
β is a reference angle of wave obliquity, or the angle between the incident wave crest (assumed linear) and the control line, joining the upcoast diffraction point to a point on the near straight beach, namely the downcoast control point. The radius
R to any point on the bay periphery in static equilibrium is angled
θ from the same wave crest line radiating from the point of wave diffraction upcoast. The three C constants, generated by regression analysis to fit the peripheries of the 27 prototype and model bays, differ with reference angle
β. Their analytical expressions have been provided by [
14]:
In contrast to the log-spiral, the parabolic equation origin coincides with the diffraction point and therefore the equation is directly related to wave direction. However, it is affected by a drawback: the uncertainty of locating the downdrift control point. Despite that there are several interpretations of the downcoast control point of the parabolic bay shape equation [
14], it is a considerable limitation which inhibits the application of the parabolic model in designing new beaches.
The hyperbolic-tangent shape model was derived by [
13] through the analysis of 46 beaches around Spain and North America. Its equation is defined in a relative Cartesian coordinate system (
Figure 3) as:
where:
= Distance across shore [m]
= Empirically determined dimensional coefficient [m]
= Empirically determined dimensional coefficient [1/m]
= Empirically determined coefficient [dimensionless]
= Distance alongshore [m]
The x-axis is parallel to the general trend of the shoreline while the y-axis points shoreward; the origin of the coordinates is placed at the point where the local tangent to the beach is perpendicular to the general trend of the shoreline. The hyperbolic tangent curve is symmetric with the x-axis and produces two asymptotes, found at
. The line
indicates the location on the shoreline which is no longer under the influence of the headland. The parameter
controls the magnitude of the asymptote (distance between the relative origin of coordinates and the location of the straight shoreline),
is a scaling factor controlling the approach to the asymptotic limit and
controls the curvature of the shape. These unknowns were found by using trial and error and an optimisation procedure that minimises rms errors. The authors of [
13] found the following relationships:
.
The model indirectly considers wave direction and, initially, it does not correlate the diffraction point with the origin of the hyperbolic tangent.
In 2018, [
14] presented a development of the hyperbolic-tangent equation, establishing a relationship between the existing hyperbolic-tangent shape equation with the wave diffraction point. The authors used a database of case studies comprised of 46 beaches in Spain, Southern France and North-Africa, determining the coefficients
a,
b and
m. The value of
m is about 0.496, supporting the original findings of [
13]. Moreover, the authors estimated a correlation between
a and
b, and the hyperbolic tangent shape equation became:
Overall, the research’s relevance has been constituted by two relationships which link the hyperbolic tangent profile and location of diffraction point. This makes the model easy to use, both in order to verify the bay’s stability and design a new beach. In the former case, starting from the distance between the location of diffraction point and asymptote, it is instantaneous to locate the origin of the equation and obtain the hyperbolic tangent profile; in the latter case, based on the shoreline advancement established in the illuminated zone, the location of the new diffraction point (i.e., new headland’s tip) can be directly carried out. The relationships, which correlate parameters
c,
d and
a shown in
Figure 4, are:
It is worth noting that the hyperbolic-tangent model could be applied not only to bay in static equilibrium conditions, but also to fit a bay under non-equilibrium conditions [
3]. Among the three aforementioned models, the parabolic model prevailed over the others since it was the only one that used the wave diffraction point as the origin of the co-ordinate system, ensuring that the effect of relocating the diffraction point can be assessed [
14]. Despite that, in comparing the origin of the three coordinates system in studying a headland bay (using a computer program based on a trial and error approach), [
25] found that the parabolic equation’s origin was located in the middle of the ocean (
Figure 5). This outcome weakens the strong point of parabolic model, revealing an uncertainty about its alleged robustness.
On the whole, regardless of models’ possible limitations, their description here presented shows how none of them model the relationship between wave characteristics and beach shape, but instead simply characterize the geometry of the headland bay profile. Nevertheless, over the last few decades, further investigations have been carried out towards the comprehension of the correlation between crenulate-shaped bays and wave characteristics. Particularly, several works focused on wave directional spreading, demonstrating that the crenulate-beach plan-form suffers under its influence. The authors of [
26], searching for a new criterion to locate the downdrift control point of the parabolic shape equation, observed that the broader the directional spreading the farther the position of the downdrift control point. This behaviour indicates that, with a wider directional spreading of waves, a greater area could be influenced by wave diffraction behind the coastal structure (i.e., the transition zone increases).
During their numerical experimentation, [
27] observed that a narrow directional spreading produces a more curved and asymmetric planform, confirming that the crenulated-beach behaviour is related to the modal directional spreading of incident waves [
26,
28]. Furthermore, [
28] evaluated how the degree of directional spreading influences diffraction’s effect on the headland shape, observing that diffraction appears less important with a high directional spreading degree. The latter outcome is in accordance with [
29], which demonstrated that diffraction shapes bay morphology only for wave climates with restricted approach angles.
Although important progress has been made, the effective correlation between wave forcing (i.e., diffraction and refraction) and static equilibrium profile has not been fully examined and understood. In this work we attempt to understand this issue, starting from the commonly engineering assumption that an equilibrium profile follows the wave front trend, with the aim of verifying the validity of this statement.
3. Numerical Experiments
The primary aim of this research is to investigate the relationship between the headland-bay static equilibrium profile and wave propagation characteristics; to do so, numerical modelling of a single headland case study has been implemented. Moreover, the possible influence of wave direction and wave period has been verified on this researched correlation. Several experiments have been carried out using the wave driver BW of MIKE 21 (Danish Hydraulic Institute) [
21], which is based on the numerical solution of the time domain formulations of Boussinesq-type equations [
30,
31,
32,
33]. The BW module is able to reproduce the combined effects of all important wave phenomena, among them diffraction and refraction, which play an important role in headland-bay equilibrium profile shaping.
3.1. BW Model of MIKE 21 (DHI)
The Boussinesq model is capable of reproducing the combined effects of important wave phenomena, such as diffraction and refraction. Therefore, it has been used in numerical experimentation to obtain diffracted wave fronts generated by headland.
MIKE 21 BW includes two modules, the 2DH wave model and 1DH wave model, both based on the numerical solution of time domain formulations of Boussinesq-type equations, which are solved using a flux-formulation with improved linear dispersion characteristics. The enhanced Boussinesq equations were originally derived by [
30,
31], making the modules suitable for simulation of the propagation of directional wave trains travelling from deep to shallow water. Moreover, it contains wave breaking and moving shorelines, as described in [
32,
33,
34]. The 2DH BW model has been used in this work.
The enhanced Boussinesq equations are expressed in terms of the free surface elevation,
ξ, and the depth-integrated velocity-components,
P and
Q. The continuity and the momentum-conservation equations read:
with the dispersive Boussinesq terms
and
, defined by:
where
d is the still water depth;
g is the gravitational acceleration;
n is the porosity;
C is the Chezy resistence number;
α is the resistance coefficient for laminar flow in porous media;
β is the resistance coefficient for turbulent flow in porous media;
B is the Boussinesq dispersion coefficient; the terms
Rxx,
Rxy and
Ryy indicate the incorporation of the wave breaking by means of the surface roller model [
35].
In MIKE 21 BW, the waves may either be specified along open boundaries or be generated internally within the model through the generation line; the latter must be placed in front of a sponge layer absorbing all outgoing waves. Moreover, porosity (e.g., to model partial transmission through porous structures) and sponge layers can be used on an ad hoc basis. At open boundaries, either a level boundary, namely wave energy given as time series of surface elevation, or flux boundary, where flux density is perpendicular to the boundary, can be set.
3.2. Methodology
The method followed involves comparing the diffracted wave fronts generated by the BW simulations with static equilibrium profiles sketched out by the hyperbolic-tangent model [
13,
14]. The procedure obeys to the following steps:
generating the diffracted wave front through the BW module;
extrapolating the wave front from the model;
starting with the illuminated zone of the front, drawing out the hyperbolic-tangent profile.
The third step is directly related to one of the major findings in the field of static equilibrium bays, namely that equilibrium crenulated beaches tend to align transverse to the direction of dominant waves. The main assumptions in applying the methodology here proposed are that (1) the equilibrium headland-bay planform follows the wave front; and (2) that the predominant wave direction is perpendicular to the straight area of the bay. In this way, in order to obtain the hyperbolic tangent plan-form of every model configuration, we supposed that the asymptote of the hyperbolic tangent matches the illuminated zone of the wave front, out of the shadow zone where the influence of diffraction is negligible. Therefore, once the distance between asymptote and diffraction point is measured (distance
c in
Figure 5), the origin of the hyperbolic tangent is automatically obtained through Equations (8) and (9), and the x,y coordinates of hyperbolic tangent profile are achieved through Equation (7) (
Figure 5). It is worth pointing out that, in order to compare equilibrium static profile and wave front, the latter must be able to expand without any kind of physical interference (e.g., presence of additional obstacles).
Additionally, in order to investigate the influence of wave characteristics on the relationship between wave fronts and equilibrium profiles, we performed different numerical simulation scenarios by varying wave direction, wave period and refraction conditions. For each scenario, more wave fronts were extrapolated from the BW module, each progressively further away from the headland tip, in order to examine the influence of dimensionless distance (c/L) on the researched correlation.
Finally, before describing the model set up, a clarification regarding the choice of the headland bay shape model is necessary. The two possible models to be implemented to sketch out the static equilibrium profile were the hyperbolic-tangent model [
13,
14] and the parabolic model [
15,
16,
17,
18] (the logarithmic spiral model [
12] has been rejected a priori given its difficulty in practical application). Thus, the same methodology has been implemented to both the hyperbolic model and parabolic model. The results proved that there was no difference: there is no change, regardless of the model adopted. Nevertheless, given the uncertainty due to the determination of the downdrift control point of the parabolic profile, the hyperbolic tangent model [
13,
14] turned out to be the best static equilibrium model to be compared to the BW simulations wave fronts.
3.3. Boussinesq Wave Module Set Up
The primary aim of this research is to investigate the relationship between headland-bay static equilibrium profile and wave propagation characteristics (diffraction and refraction). Therefore, we created different numerical simulations scenarios, in order to generate different wave fronts, by varying wave direction, wave period and refraction conditions.
First, to evaluate refraction conditions, two bathymetry configurations were been investigated. The first one (called “
gentle slope”) is characterised by a linearly varying cylindrical bathymetry, with a gentle bottom slope of 1/100 and a diffraction point, modelled through a breakwater, located at a water depth of 3 m. The water depth at the offshore boundary is 20 m (
Figure 6a). A gentle slope has been adopted to allow an expansion of the wave fronts without the influence of bottom abrupt raising. The second configuration (called “
flat bottom”) exhibits a bottom with a constant water depth equal to 10 m (
Figure 6b).
Secondly, numerical experiments have been conducted considering the possible effect of wave period and wave direction; conversely, the influence of wave height was not explored. This is because the correlation function describing the relationship between static equilibrium profiles and wave characteristics, at this first stage of the research, is obtained through a “linear approach” based on simple geometric consideration. Therefore, the only parameters that could geometrically affect the aforementioned correlation are period, wave direction and bottom inhomogeneity. Conversely, wave height does not represent a variable of the problem. Three different wave directions have been investigated; one is normal to the breakwater, while the other two are angled of 25° and −35° with respect to being perpendicular to the structure. For each direction, three wave periods have been tested: 5.9 s, a typical wave period of the Mediterranean wave climate; and 10 s and 15 s, which simulate swell conditions. The wave height used was 0.8 m for each numerical scenario. Since the predominant wave shapes the crenulated beach, for the numerical experimentation the value of wave height has been chosen equal to the LDR equivalent wave height value for the case study of Bagnoli bay (see
Section 5.2). Simulations have been implemented using regular waves. It is important to underline that the effect of breaking waves was not considered in the trials. For each bathymetry configuration (“gentle slope” and “flat bottom”), the scenarios analysed are summarized in
Table 1.
Models are made up of a fine grid (square cells with grid spacing 3 m), upon which orientation coincides with wave direction; the time step used is 0.1 s. The wave generation line has been used and wave absorbing sponge layers have been applied at the model boundaries. Near the land, one sponge layer has been applied in order to avoid the occurance of wave reflection which could influence the expansion of the wave fronts. Geometric characteristics are summarized in
Table 2, and are valid for both bottom configurations.
6. Numerical Modelling of the Bay and Results
Numerical modelling of the bay has been performed through the wave driver MIKE 21 BW (DHI), which allows for obtaining the growth of the wave fronts leeward of the headland. In order to follow the method used in the numerical trial, the bay has been modelled in such a way that wave fronts could extend leeward the headland without any interference: the real bathymetry of the area has been employed until the diffraction point location, where the water depth is about 3 m, then a gentle slope of 1:100 has been adopted; the coastline has been shifted and cut landward and a sponge layer has been used, in order to avoid reflection phenomenon (
Figure 17).
Numerical simulations have been performed wave detecting the equivalent direction of 205° N, since it is representative of the effects of the entire wave climate on long-shore sediment transport [
41]. As concerns wave period, the average measured omni-directional peak period (5.9 s) has been investigated.
Simulations have been carried out using regular waves, which wave height has been set at 0.8m; the breaking phenomenon has been neglected. Fine grid has been used (square cells with grid spacing 3 m), for which the orientation coincides with wave direction, with a time step of 0.1. A wave generation line has been used and wave absorbing sponge layers have been applied at the lateral boundaries. Grid geometric characteristics are summarized in
Table 3.
Comparison between equilibrium profile and wave fronts carried out through BW model, shown in
Figure 18b, demonstrates that wave fronts’ direction in the illuminated zone matches shoreline orientation, confirming the goodness of the LDR equivalent wave to represent the entire wave climate in the act of shape the beach planform. Specifically, simulations results confirmed the behaviour observed in the previous section: the shoreline planform of the headland bay is placed on more wave fronts, the downdrift section overlaps the wave front while the curved section is placed on the wave front closer to the headland (
Figure 18a).
In order to attain an effective validation of the
wave-front-bay-shape equation, Equation (15) has been applied. The minimum distance between the diffraction point and the asymptote has been measured (distance c in
Figure 19a) and it has been standardised with respect to the local wave length,
c/L. Applying Equation (15), we obtained the value of the shift,
s, which represents the distance needed to superimpose the whole wave front on the shoreline. After shifting the wave front, as shown in
Figure 19b, it is perfectly superimposed on the bay shoreline, thus verifying the correlation found between wave fronts and shoreline profile, reached by Equation (15).
7. Conclusions
The main purpose of present paper is to establish a relationship between static equilibrium profiles of crenulated bays and wave front shapes and, more generally, understanding the correlation between bay shoreline response and wave forcings. The peculiar equilibrium plan-form (static or dynamic) of headland bays is assumed at a relatively long-term scale (e.g., annual to decadal) as a response to predominant wave direction; the downdrift segment of the bay, in fact, tends to align to the dominant direction of incoming waves, while, the curved up-coast segment is modelled by the diffracted wave fronts. In this regard, it is commonly believed by engineers that equilibrium beach profiles follow the wave front trend; however, this has been never fully proved so far. The most important works concerning headland bay beaches provided by the literature are empirical models, which simply describe the geometry of static equilibrium bays shape, neglecting the acting physical processes that sculpture the coast. In fact, they are merely of an empirical kind, lacking in a further insight on relationships between incident wave characteristics and beach shape. Therefore, in the lack of a model in which plan-shape is strongly correlated with wave characteristics, the project design of a new headland bay beach (necessary if the headland control practice has to be implemented) could, as a result, be extremely challenging.
To assess a possible correlation, numerical experiments have been carried out using the MIKE 21 Boussinesq Wave Module (BW), where wave fronts have been compared to the hyperbolic-tangent equilibrium profile, analysing the influence of wave direction, wave period and refraction phenomenon. Results proved that equilibrium profiles are located seaward compared to their relative wave fronts. Hence, designing a new beach following a wave front trends ensures a shoreline accretion in the sheltered zone (i.e., leeward to the headland). Additionally, a correspondence function, called the “wave-front-bay-shape equation” has been established, offering an easy application to engineering uses due to the simple geometric interpretation of its controlling parameters. The function seems to indicate that equilibrium beach profiles of a single headland bay correspond to a simple translation of wave front normal to the propagation line at the headland tip.
Moreover, the application of the “wave-front-bay-shape equation” to the case-study bay of the Bagnoli coast (south-west of Italy) has been performed. The numerical model has been set up, and the LDR equivalent wave concept has been used to embody the dominant wave attack that rules the long-term evolution of the little bay. Results confirm the behaviour observed from numerical experiments outcomes; the “wave-front-bay-shape equation” has been successfully verified, thus confirming that a correlation between equilibrium plan form and wave fronts can be found.
Nevertheless, the research still stands at a primary stage, and requires improvement and accurate verification in future research works in order to develop an enhanced guidance which could help in engineering and morphological practice.