Sediment accumulation in embayments controlled by bathymetric slope and wave energy: Implications for beach formation and persistence

High energy, rocky coastlines often feature sandy beaches within headland‐bound embayments. Not all such embayments have beaches however, and beaches in embayments can be removed by storms and may subsequently reform. What dictates the presence or absence of an embayed beach and its resilience to storms? In this paper, we explore the effect of offshore slope and wind conditions on nearshore sediment transport within idealised embayments to give insight into nearshore sediment supplies. We use numerical simulations to show that sand can accumulate near shore if the offshore slope is >0.025 m/m, but only under persistent calm conditions. Our modelling also suggests that if sediment in an embayment with an offshore gradient steeper than 0.025 m/m is removed during a period of persistent stormy conditions, it will be unlikely to return in sub‐decadal timescales. In contrast, sediment located in embayments with shallower gradients can reform swiftly in both calm and stormy conditions. Our findings have wide implications for contemporary coastal engineering in the face of future global climate change, but also for Quaternary environmental reconstruction. Our simple method to predict beach stability based on slope can be used to interpret differing responses of embayments to periods of changing coastal storminess such as the medieval climate anomaly‐little ice age (MCA‐LIA) transition. © 2018 The Authors. Earth Surface Processes and Landforms published by John Wiley & Sons Ltd.


Introduction
Numerous factors control the form and persistence of sandy beaches, including sediment supply, incident wave direction, alongshore and cross-shore sediment transport mechanisms, and tidal range (Masselink and Short, 1993;Masselink and Pattiaratchi, 2001;van Rijn et al., 2003;Aagaard et al., 2004;Cooper et al., 2004;Austin et al., 2009). A number of authors have explored the morphological response of existing embayed sandy beaches to wave and storm attack (Masselink and Short, 1993;Short, 1996;Forbes et al., 2004;Jackson et al., 2005), and embayment morphology is known to be sensitive to wave driven alongshore (Ratliff and Murray, 2014;Hurst et al., 2015) and cross-shore (Harley et al., 2011) sediment fluxes. The patchy recovery of embayed beaches in the face of storm action has been also noted recently in several studies (Loureiro et al., 2009(Loureiro et al., , 2016Roelvink et al., 2009). The long-term presence or absence of beaches has received far less attention, however. Recently landscape-scale models of beach and headland erosion have focused on sediment supply from headland erosion to maintain beaches (Limber and Murray, 2011), but less attention has been paid to the sediment transport dynamics that may promote or inhibit beach formation and stability.
Beaches have important economic and cultural value to coastal communities. Not only is the desirability of living by the coast a major economic factor, with associated tourism and recreational activities, but also beaches and coastlines hold many cultural and archaeological sites. Understanding what drives their formation, stability and distribution is crucial for effective management strategies and adaptation in the face of future climate change (Zhang et al., 2004;Dawson, 2013). Beach sediment is eroded during storm events and often lost offshore, but beaches can recover these sediments gradually over seasonal to decadal timescales (Harley et al., 2015;Scott et al., 2016).
Beach formation is dependent on the availability of sediment, as well as its movement and residence time at the nearshore. Sediment can be derived from three basic sources: onshore (e.g. shoreline-derived sediment from coastal erosion or delivered by rivers), offshore (e.g. glacial or fluvio-glacial sediments, other sub-marine sediment and biogenic materials, re-worked by rising sea levels) and alongshore (i.e. supplied by alongshore drift). At specific locations, sediments from these different sources are reworked by the interplay of currents, waves and wind with topography (Castelle and Coco, 2012;Maspataud et al., 2009).
Offshore sediment supply to replenish beaches is transported nearshore and is often deposited as sediment berms and sand bars, the formation and stability of which can be impacted by storm events (Ruessink et al., 2016). High-energy coastlines without headlands lack barriers to the lateral movement of sediment, thus beaches can be formed or replenished wherever alongshore drift is active and there is sufficient sediment supply. In contrast, the irregularities of high-energy, headland-dominated coastlines inhibit alongshore sediment flux (Short, 1996;Limber and Murray, 2011;da Silva et al., 2016), allowing beach material to accumulate between headlands to form embayed beaches (Limber and Murray, 2011). Aeolian erosion and deposition, particularly in these storm-prone high-energy environments, can be a driver of beach destruction and formation without oceanic interaction. For example, some beaches in sheltered embayments may be due to interactions with dune systems, with cycling of sediment replenishing beaches from dune material, and vice versa (O'Connor et al., 2007), but in many cases dunes are absent and a beach-hinterland cycling model cannot explain phases of beach formation and destruction.
Despite similarities in sediment supply, exposure to wave and tides and other forces as described above, we observed that not all embayments along a coastline will contain beaches ( Figure 1).
Our aim is to better understand what controls the input of offshore sediment to beaches on high energy, headland dominated coasts and what controls the formation, stability and replenishment of these beaches. In this paper we combine morphological analysis and coupled morphodynamic and hydrodynamic exploratory modelling to investigate how offshore slope and wind regime, two driving factors behind coastal sediment transport, can affect migration and accumulation of sediment in the nearshore environment.

Approaches and methods
We adopted a research strategy that combines empirical observation and numerical modelling. First, we analysed bathymetric data from examples of headland-dominated, high-energy coastlines to highlight an empirical relationship between local offshore slope and the presence or absence of sandy beaches in headland-bounded embayments. We sought to understand this relationship through numerical modelling of embayments with different offshore slopes and wind conditions.

Observations of offshore gradient and the presence of beaches
We used bathymetric charts compiled by the United Kingdom Hydrographic Office (UKHO) to quantify offshore slope in a range of locations ( Figure 2). Formerly glaciated coastlines with high-energy, headland-dominated coastlines were chosen for analysis, because extensive glaciation leaves a legacy of unconsolidated offshore sediment that could form a ready supply of material for beach formation (Clark et al., 2012). Bed elevations were determined at 1 km from shore roughly perpendicular to the orientation of the shoreline in the centre of the embayment. These were used to calculate the average offshore slope. We fully acknowledge the seemingly arbitrary nature of the 1 km distance, however, this distance was chosen to avoid a bias towards small-scale nearshore features, such as shore platforms, that are poorly resolved on most Admiralty chart data.
At each site, we selected stretches of headland-dominated coastline and measured average offshore slope in every embayment within that stretch. Coastlines were initially identified by availability of high resolution Admiralty charts (greater than 1:75 000). Sites were not filtered by onshore factors. We removed embayments that had skerries or barrier islands within 1 km of the shore. Embayments with major inputs of sediment from terrestrial sources (such as those containing river deltas) were also omitted, to ensure that sediment supplies moving across the offshore slope was the principle factor under consideration (see Supporting Information).
Google Earth imagery was used to determine the presence or absence of a sandy beach. A qualitative analysis of the colour and texture of beaches on aerial imagery supported by any other available photographic evidence (for example, Panoramio photographs in Google Earth, and Google Image searches for specific embayments) was used to determine if beaches were sandy or not. A total of 239 embayments were analysed, of which 119 contained sandy beaches. These data were grouped based on 0.05 m/m increments in offshore slope derived from admiralty charts. We calculate the probability of finding a beach as the number of beaches in the slope bin divided by the number of sites analysed in that bin. We find that whereas the majority of sites with slopes less than 0.025 m/m had beaches, far fewer beaches were found at sites with slopes greater than 0.025 m/m. For slopes greater than 0.025 m/m, the probability of finding a sandy beach decreased with increasing bathymetric slope ( Figure 3). During sampling, few suitable embayments were identified with gradients above 0.034 m/m and below 0.009 m/m, thus these end members bins should be treated with caution.
We acknowledge that compounding factors exist that may impact the relationship presented in Figure 3. Specific morphological factors such as sediment composition and headland amplitude, as well as hydrodynamic factors such as local wave and tide conditions can all impact the accumulation of sediment in embayments suitable for beach formation. However, this association was of sufficient clarity for us to undertaken exploratory morphodynamics modelling to better understand the controls on this relationship.

Numerical modelling
To better understand our semi-quantitative observations and explore the potential role of bathymetry and coastal sediment transport in dictating the presence or absence of embayed beaches, we performed simulations of the evolution of idealised headland embayments using the MIKE21 coupled hydrodynamic and morphodynamic model. This model was chosen because more than 20 years of development has produced a robust and reliable tool for coastal engineering studies (Siegle et al., 2004;Manson, 2012;Houser, 2013;Vicinanza et al., 2013;DHI, 2014;Contestabile et al., 2015). The two-dimensional hydrodynamic (HD) and sand transport (ST) modules were coupled, and hydrodynamics were forced by the spectral wave (SW) module, which was run independently to calculate wind-generated waves and wave radiation stresses as an input for the HD module. The SW includes radiation stresses which ensures conservation of momentum for waves breaking nearshore; this module is necessary to compute nearshore currents and thus nearshore sediment transport (DHI, 2014). Calculated radiation stresses for model scenarios are added as inputs into the HD module.
The MIKE21 ST module moves sediment in the domain by calculating the shear stress generated by water motion, using sediment transport theory for both bedload and suspended load developed by Engelund and Fredsøe (1976) and Fredsøe et al. (1985). The module treats sediment as the total load (i.e. bedload + suspended load). The water motion leading to shear stress and sediment transport combines wind-generated waves, tides and currents, which are calculated by the HD and SW modules. MIKE21 calculates shear stress by combined wave and tidal currents, which cannot be determined separately. However, we present here shear stress calculations for waves, as this is the dominant force acting on sediment motion above the wave base in the model. Sediment motion is initiated when the shear stress τ exceeds the critical shear stress τ cr required for incipient particle motion.
Dimensional shear stress τ is defined: where ρ ω is fluid density, C f is the wave friction factor and U b is the horizontal (i.e. cross-shore) wave orbital velocity at the bed. Soulsby and Clarke (2005) calculated the wave friction factor by: where T is the wave period, and D 50 is the median grain size. The HD module does not directly output the shear stresses or wave friction factors, although these can be calculated per time step from the horizontal wave orbital velocity outputs following these equations (Soulsby and Clarke, 2005;DHI, 2014). The critical shear stress τ cr is dependent on particle size and density, and is determined by: Where θ cr is the non-dimensional Shields parameter, ρ s is sediment density and g is the gravitational constant (Shields, 1936;Madsen and Grant, 1976;Paphitis, 2001). Where τ > τ cr in the domain, the ST module computes sediment transport following Engelund and Fredsøe (1976) for bedload and Fredsøe et al. (1985) for suspended load. Both bedload and suspended load transport are solved deterministically. Timeaveraged bedload transport is solved based on current outputs from the HD module. Suspended load is continuously active within the model domain provided that sediment fall velocity (as a function of suspended sediment concentration) is exceeded by current velocities computed by the HD module. The model sums bedload and suspended load transport to calculate the total sediment load and updates the bed morphology accordingly each time step.

Model domain
We designed an idealised embayment formed by a rectilinear grid, 2 km in length and 1 km in width, created using the MIKE21 domain generator. We chose these dimensions to be as representative of embayments found along high-energy coastlines from our Admiralty chart analysis ( Figure 4).
We set a non-erodible land boundary for the shoreline, and two land boundaries representing non-erodible headlands for the lateral edges of the domain (see Figure 5). The open boundary was set so that the spatial change in sediment flux across the boundary was zero. This allowed sediment to both enter and leave the domain as demanded by the changing hydrodynamic conditions, without suddenly depositing or eroding material at the boundary, avoiding a glass wall effect (Keen et al., 2003;Manson, 2012). We chose this boundary condition to avoid artificial accumulation of sediment against the boundary. Table I lists the model initial conditions. Note that these boundary conditions are specific to embayments that do not have a significant terrestrial sediment input.
Non-erodible bathymetry (the hard bedrock substrate) was created as a uniformly sloping plane where the seaward side (1 km from shore) was set to depths ranging from 20 m to 30 m, to generate shoreface slopes ranging from 0.02-0.03 m/ m. We added low-amplitude noise to the nearshore bathymetry with an average amplitude change of 0.15 m ( Figure 5). This is intended to reflect minor morphological heterogeneity to the underlying bedrock surface by simulating small irregular undulations in nearshore contours. These perturbations are 'smoothed out' after initial model time steps and final model results are not sensitive to the amplitude of the noise imposed.

Initial conditions
We added a uniform 1 m thick sand layer across the entire domain, representing an initial sediment volume in the domain of 2×10 6 m 3 . The properties of this sediment layer were generated in MIKE21 using the Q3D sediment generation table. This utility produced a file containing sediment properties that were used by the ST module and applied for each time step. Table II lists the sediment properties chosen for the modelling exercise. There is no simulation of the actual presence of a beach above the water line as MIKE21 cannot resolve mesh elements that become 'dry' (i.e. at low tide where a beach may form), but sand bars form just below minimum water level, which acts as a ready sediment supply for beach building. Our aim here was to explore climate and slope effects on the shoreward migration of coastal sediment.   We chose a meso-tidal range of 2.56 m for the modelling, derived from the tidal range as recorded in Lerwick, Shetland, UK. This is a typical high-energy headland-dominated coastline with a long tidal gauge record. Model sensitivity analysis shows the model to be insensitive to tidal range (see Supplementary information).
To remove bias of specific storm surges or other tidal swells, we based the tidal record on the Lerwick data. Rather than using gauge records, this was generated using tide-prediction software JTides (JTides, 2017) (To capture the impacts of wave action on sediment transport, a time step of 30 min was used for the tide. A tidal record of 1 year was generated (the maximum the program allows), and extended to create a 10 year record. We chose a simulation time of 10 years to provide a sufficiently 'long-term' evolution of an embayment, capturing the influence of multiple years of storm impacts on the coastline. This allows an understanding of beach formation, persistence and recovery framed by a timescale on which economic uses of the coastline rely, whether modern or historical.
The specified significant wave height (SWH) and the peak wave period (T) are modified by wind-generated waves in the HD and SW module, and act as initial conditions for the wave field. Wind fields specified in the HD module modify surface waves as they transform across the bathymetry, modifying SWH and period each time step. SWH and T values were derived from average values over a year (2011) from the closest sea buoy to Shetland (NOAA, 2016) (see Table II). While storm surges are important for subsequent sediment transport, we do not include any specific surges beyond that which is calculated as part of the SW results as the purpose of the experiment is to understand nearshore sediment transport on a long-term basis, rather than any specific storm.
Radiation stresses produced by wave action are governed by the output of the SW module and applied to the HD module for the simulation of sediment transport. Radiation stresses are calculated based on SWH, which is subsequently modified each time step due to wind-sea and swell conditions. Thus there is no specific swell-wave significant wave height chosen. Peak spectral wave period (the wave period corresponding to the maximum wave energy level in the wave spectra) is an output of the spectral wave module based on an initial peak wave period. Wave direction at the offshore boundary is set as perpendicular to the model shoreline (0°) to represent fetch-unlimited waves impacting the shoreline without interference from refraction around the headlands (or islands). Wind-sea and swell conditions are calculated in the spectral wave module, with the JONSWAP formula added. While this is fetch-limited by definition, the maximum fetch length specified in the model is 1000 km, thus the formula calculates an essentially 'unlimited' fetch for our purposes. Wind-sea and swell conditions are derived from the SW module and applied to the HD module.
We split wind forcing into two categories, calm and stormy. Median wind speeds on Shetland (as stated, a typical highenergy coastline prone to storminess) are~7.36 m/s (average monthly mean 1930-2010), and so conditions were chosen to reflect this. Thus calm conditions were specified to range from 1 to 15 m/s, and stormy conditions to range from 1 to 60 m/s. The maximum value of 60 m/s (active for only 6 hours out of 10 years of model simulation) was chosen to represent persistently stormy conditions on the coastline as it is the median of the highest wind gusts recorded in Shetland in each of the past 30 years, which range from 45 m/s to 77 m/s (Shetland Islands Council, 2011).
We set wind direction to be randomly selected between north west (325°) and north east (45°) every 6 h. Figure 6(a) illustrates the running mean wind speed distribution (over a model 3 day period) applied in the model over a typical model 'year' in 6-hourly increments, for calm and stormy scenarios. While this method does not take into account seasonal storminess, we chose this to subject a headland embayment to constant storminess and thus determine the physical effect these persistent conditions have on sediment transport onshore, regardless of storm season. Onshore winds (i.e. those blowing from onshore to offshore) were ignored, as were across-shore winds (i.e. winds parallel to shore). This is a limitation, but we deliberately built this into the experiment in order to isolate the effect that purely offshore wind-generated waves and subsequent wave energy has on nearshore sediment transport (Van Donk et al., 2005). Wind speed randomly varied at the same temporal frequency based on a Weibull distribution of the bounds listed above (Figure 6(b)), common of wind speed distributions at many coasts (Tuller and Brett, 1984;Van Donk et al., 2005;Kidmo et al., 2015).

Sensitivity analyses
We carried out the sensitivity analyses on the model to lend confidence that our results are not unduly influenced by choice of model parameters (Table III). Sensitivity analyses can be found in supporting information.
While both bedload and suspended load are operating in the model and reported as total load, both Bagnold (1966) and Bowen (1980) demonstrated bedload transport was the dominant process nearshore for sand particles of 250 μm in diameter when horizontal orbital velocities are less than 0.35 m/s. Mean horizontal orbital velocity in calm scenarios for all slopes is less than 0.1 m/s, and less than 0.3 m/s in stormy scenarios for all slopes. MIKE21 calculates model results based on mesh elements in the domain, essentially model 'cells', which can then be interrogated once the simulation is complete. Both numerical and graphical representations of mesh elements are output. These were interpolated from the model results mesh and exported to a 10 m resolution grid in order to perform differencing in bed elevation and to visualise results.

Modelling results
We explored the morphological evolution of idealised embayments with contrasting bathymetry and differing wave forcing conditions. Sand bars form just below mean low water for both shallow and steep slopes under calm conditions, although on the steeper slope sediment deposition is patchy, with thinner sediment deposits towards the centre of the model domain. MIKE21 does not resolve sediment accumulation within the swash zone where model cells become 'dry' as the tide ebbs (approximately 150-200 m from the shore, depending on offshore slope), but sediment for beach formation is available offshore if the bar is present. Under both calm and stormy conditions, the sand bar formed after one month of model simulation time and stayed relatively uniform throughout the 10-year simulation time.

Planform sediment patterns
A marked difference in shallow and steep slopes was observed in terms of sediment distribution forced by stormy  conditions. On shallow slopes, sand bars form in both calm and stormy conditions but sediment is thinner in sand bars formed in stormy conditions. On steep slopes, sand bars also form in calm conditions however, no sand bar formed under stormy conditions, with sediment instead rotated towards the headland walls and offshore into deeper water. Figure 8 shows swath profiles of sediment thickness at the end of the model runs. In both calm and stormy scenarios, sediment has aggraded toward shore, whereas sediment has been lost (i.e. is less than the initial 1 m thickness) in the seaward side of the model domain. The steep drop-off of sediment thickness close to shore is representative of the tidal range in the model, and thus sediment transport cannot be resolved in this part of the domain as the cells go 'dry' as the tide ebbs. We find maximum sediment thickness at offshore distances just below mean low water (MLW). In both calm and stormy scenarios, there is a distance from shore where sediment thickness increases from the initial 1 m, varying between 450 and 550 m offshore dependent on slope. These locations are roughly coincident with wave base, 14.4 m in our model set up, defined as half the wavelength (Peters and Loss, 2012), varying slightly between stormy and calm conditions. We therefore restrict our analysis of volumetric sediment flux to the shoreface landward of this depth where there is significant change in bed elevation and sediment thickness. Figure 9 shows volume change over time for each slope scenario, for calm and stormy conditions in the shoreward half of the domain. The sharp drop in volume seen between initial conditions and the first model time step in both calm and stormy scenarios is caused by initial sediment mobilisation in the model domain.

Volume flux
Under calm conditions, sediment volume generally increases over time in the shoreward half of the domain (Figure 9(a)), consistent with the planform sediment accumulation seen in Figure 8. The increase in sediment accumulation is fed by the zero sediment flux boundary in the model. Figure 10 shows the mean sediment flux within the model domain, averaged per time step, as a function of bathymetric slope. Under calm conditions the rate at which sediment accumulates in the model domain decreases with steeper bathymetric slope, particularly when the slope exceeds 0.023 m/m.
Conversely, under stormy conditions, sediment volume in the model domain decreased through time (Figure 9(b)) and the rate of sediment loss was insensitive to the bathymetric water. Sediment depletion rates are relatively insensitive as slope increases under stormy scenarios. An interesting impact can be seen in Figure 9(b) just before year 9 of the stormy simulation, in which volume flux becomes positive for slopes of less than 0.026 m/m due to wind speeds falling to 15 m/s for an extended period of time, but then return to negative flux when stormy conditions resume. For slopes of 0.026 m/m and steeper, no positive flux is seen. This negative sediment flux is exhibited for all modelled scenarios but a weaker relationship exists between mean sediment loss per time step and increasing slope. Rates of sediment loss per time step are relatively stable as slope increases between 0.022 m/m and 0.028 m/m ( Figure 10).
If we further explore the change in sediment flux observed in Figure 9(b) when climatic conditions briefly switch from stormy to calm, Figure 11 shows the impact on shoreward volume of switching climatic conditions halfway (i.e. 5 years) into the simulation: Change in volume flux as slope increases is sensitive in a nonlinear fashion to a stepwise change from calm to stormy conditions (Figure 11(a)). For shallower slopes (0.020 m/m to 0.024 m/m), volume flux response time is almost immediate once stormy conditions take effect, with volume flux switching from positive to slightly negative (section i), Figure 11). A threshold gradient appears to be 0.025 m/m, consistent with the Admiralty chart analysis, from which steeper slopes all share a similar volume flux response time of approximately 3 months from the onset of stormy conditions to the change between positive and negative volume flux (Figure 11(a); section (i). As the scenarios progress, slopes shallower than 0.025 m/ m begin to exhibit some volume recovery approximately 2.5 years (Figure 11(a); section (ii) from the onset of stormy conditions but at a reduced rate when compared with calm conditions. Slopes steeper than 0.025 m/m do not exhibit this recovery during the model simulation time and continue to flatline (Figure 11(a); section (iii)), suggesting no loss or recovery of sediment nearshore on this timescale. When stormy conditions give way to calm conditions (Figure 11(b)), volume flux on shallower slopes begins to increase almost immediately for slopes shallower than 0.025 m/m.

Bed elevation change
We present comparison plots of sediment accumulation as a function of the bed elevation in Figure 12 for all model time steps, with initial bed profile plotted for comparison. The standard deviation indicates the areas where there are gradients in sediment flux through time, and thus represent the areas of sediment motion in the model. Sediment becomes mobilised across the entire model domain in all cases, however, the majority of sediment motion occurs above wave base in most scenarios. Under calm conditions, significant sediment motion occurs above the wave base in all scenarios, however, during stormy conditions, some sediment motion can be seen occurring below the calculated wave base along the steeper slope, beginning to occur at the 0.025 m/m slope. This may indicate the storm wave base being deeper than under calm conditions.

Shear stress
Shear stress is not a direct output of the model, but can be calculated using Equations (1) and (2) from the horizontal orbital velocity model output. Figure 13 and Figure 14 show both absolute mean and maximum dimensional shear stress in both shoreward and seaward directions and a probability density function (PDF) for Shields values for the 0.020 m/m and 0.030 m/m slopes, for calm and stormy scenarios. The mean and max shear stresses are calculated, as for the volume flux, of the shoreward half of the domain. This is the area of significant sediment motion, above wave base, and thus is most critical to the understanding of how storminess and bed slope impacts sediment accumulation nearshore. Values reported are taken at approximately 3-day model intervals: Figure 13 and Figure 14 when viewed together suggest slope does not play a significant role in the bed shear stresses experienced in the model domains, with minimal difference throughout the model simulations between the shallow and steep slopes. Both time step plots and PDFs show the Shields values, and therefore shear stress, increasing by an order of magnitude between calm and stormy scenarios, well above the 0.05 threshold for incipient motion in the model. The PDFs overlay the Shields number in both the shoreward and seaward direction, which are seen to be roughly equal. Under calm scenarios, sediment is under entrainment <10% of model time, and >50% of model time under stormy scenarios.

Discussion
Our model results show that when an offshore sediment supply is present within 1 km of shore, the shoreward accumulation of sand-sized sediment in embayments under calm climatic conditions is relatively uniform but diminishes with increasing bathymetric slope (Figure 9(a)). With a supply of sand-sized sediment for beach material from offshore, and all other conditions equal, pocket beaches could form and persist in embayments with offshore slopes of all gradients. In contrast, under stormy conditions (Figure 9(b)) the sediment distribution becomes more variable and volume flux is negative, with net sediment loss nearshore, which over time would render a pocket beach unviable due to the starvation of nearshore sediment. If environmental conditions switch from relatively calm to relatively stormy, sediment volume change on offshore slopes with a gradient of less than 0.025 m/m begins to turn positive after approximately 2 years of stormy weather, while volume flux on steeper slopes remain slightly negative. Switching from stormy to calm climatic conditions, volume change on shallow slopes turns positive almost immediately and volume begins to accumulate nearshore, whereas steeper slopes exhibit a lag time of several months before volume change turns positive. Bathymetric slope is a major control on the nearshore sediment budget under calm climatic conditions, and thus critical to the recovery and persistence of a sandy beach after storm-induced erosion. However, stormy conditions have a stronger impact on sediment loss nearshore than slope. When climatic conditions are variable, slope controls how quickly nearshore sediment budgets recover. The accumulation of nearshore sediment berms that the model simulates are broadly consistent with recent empirical studies of beach evolution undertaken on specific embayments (Loureira et al., 2009(Loureira et al., , 2016Harley et al., 2015), as well as agreeing with flume experiments undertaken by Ruessink et al. (2016) in terms of sand bar accumulation nearshore. Our model results are also consistent with recent examples of storms that have removed beaches that have received attention in the news, including Porthleven beach in Cornwall UK in 2015 and Collaroy beach in Sydney, Australia in 2016.These beaches exist on very shallow offshore slopes, and thus our model results suggest they are liable to recover quickly once calm conditions return. Our modelling results suggest that a beach in an embayment with a steeper nearshore gradient, depleted by successive storm action, is less likely to recover its beach volume during inter-storm calm periods than a lower gradient embayment. Figure 15 is a conceptual summary of model behaviour, consistent with observed patterns of beach location from the Admiralty chart analysis.

Processes of beach formation/removal
Model results suggest that a relatively well understood concept of coastal morphodynamics, the wave base (the maximum depths at which wave orbital motion is considered significant), can be used to explain the mechanisms observed in the model results and the Admiralty chart analysis, in that the formation and persistence of sand bars, and thus beach material, is primarily dependent on wave energy.  Figure 8 shows that at most a few tens of centimetres of bed change occurs further than 450-550 m offshore compared with changes approaching a metre further onshore, depending on slope. These distances are coincident with the calculated deep water wave base (Figure 13) of 14.43 m, and change as a function of bed slope. While these figures shows minimal bed change at depths deeper than the wave base (although significant changes during initial model time steps), this does not necessarily mean minimal sediment flux. Sediment transport at depths greater than wave base are likely to be current driven, while transport above the wave base is likely to be wave dominated, consistent with the modelling results. Wave heights may be amplified due to wind forcing during stormy conditions, and thus the largest waves may agitate the bed at greater depths.

Bed shear stress
The shear stress exerted on sediment in a system, controlled by wind and wave action, governs sediment transport. Sediment is at rest until enough shear stress is imparted upon sediment to initiate motion, which is described in terms of dimensionless shear stress by Equation (3) (Shields, 1936;Paphitis, 2001;Simões, 2014). Under calm conditions, mean bed shear stress is mostly under this value, although it rises above it at multiple points throughout the model simulation ( Figure 13 and Figure 14). Under stormy conditions, mean bed shear stress is roughly an order of magnitude higher than the critical shear stress parameter for incipient bed motion, and thus is under almost constant entrainment throughout the model simulation. Patterns of planform sediment accumulation (Figure 7), show that sediment is depleted at the centre of the embayment under stormy conditions and generally rotated towards the headland as waves diffract around the headlands walls.
Sediment motion nearshore is governed in part by horizontal wave orbital velocity (Dean and Dalrymple, 1991;Chou et al., 2015;Aagaard and Hughes, 2017), in which sediment bedload will move both shoreward and seaward. Increased horizontal orbital velocity, and thus shear stress as describe above under stormy conditions carries sediment further upslope as each wave passes, but also further downslope when the orbital velocity reverses. A small volume of sediment under each wave is lost below the wave base. This volume is not replaced with subsequent waves under stormy conditions, and over time, there is a net sediment volume loss shoreward. Constant sediment entrainment under stormy scenarios explains the insensitivity of the model to an increase in tidal range. An increase in the swash zone and thus an increase in surface area along the shoreface for sediment to accumulate does not promote an increase in sediment accumulation in stormy scenarios as the sediment is removed and rotated towards deeper waters. This explanation is also consistent with the coincident maximum sediment accumulation at wave base depth (Figure 12) found in both calm and stormy conditions.

Implications
The study presented here shows an empirical relationship exists between beach formation and persistence, and average offshore slope. Our modelling results suggest that while embayments with a range of offshore slope gradients have the ability to form nearshore sand bars and thus beach material, embayments with steeper offshore slope, and those subjected to stormier conditions, are less likely to maintain a beach. These relationships have real world implications, as in periods of changing magnitude and frequency of storm events, such as from the Medieval Climate Anomaly to the Little Ice Age, or indeed from the present to future climate scenarios of the latest IPCC report (IPCC, 2013), our results highlight that there is the potential for seemingly persistent beaches to disappear and never reform. To counter coastal erosion and improve the economic and social function of coastlines, beach nourishment is often used to rebuild eroding beaches (Gopalakrishnan et al., 2011). Yet our modelling suggests that for embayed beaches where offshore gradients are steeper, long-term success is unlikely, or may require increasing volumes of beach sediment or frequency of nourishment due to significant losses of sediment offshore during storms. Where beaches are currently under threat due to increased storminess, maintenance through beach nourishment schemes may be unviable in the face of predicted climate change (Vitousek et al., 2017) .

Conclusions
Our review of bathymetric maps and satellite imagery indicated that the presence or absence of beaches on high-energy, headland-dominated coastlines might be controlled primarily by the average gradient of the shoreface. Embayments with gradients steeper than 0.025 m/m were notably less likely to have a sandy beach than those with a shallower bathymetry. We have used numerical modelling to explore the physical processes that could explain the observed relationship. Our experiments have shown that embayments can form under calm conditions on coastlines with steep offshore slopes, yet persistent stormy conditions inhibit formation of beaches in steeper sloped embayments. The threshold gradient also appears to control the recovery time of these beaches; the return to positive sediment flux on nearshore slopes (> 0.025 m/m) can take months once stormy conditions give way to calm conditions. Successive or persistent stormy conditions can completely deplete the sediment supply, particularly on steeper slopes, such that an embayed beach lost during a storm may not recover due to the loss of supply with beach material from offshore.
Under stormy conditions, shear stresses nearshore tend to be greater than the critical shear stress and sediment thus cannot accumulate nearshore to provide material for beach replenishment. This research has implications for future coastline management in terms of a new understanding of the stability of pocket beaches along headland-dominated coastlines, and the potential for this to be used as a predictive tool. The mapping of coastlines in terms of their stability can be undertaken, providing a useful resource for planners and coastal engineering efforts along socially and economically important coastlines, e.g. Fitton et al. (2016). Figure S1. Typical "calm" waves as they transform across a transect of the centre of the model domain. Figure S2. Typical "stormy" waves as they transform across a transect of the centre of the model domain. Figure S3. Wave angle change (20°and 60°) sensitivity analysis for 0.020 m/m slope, (a) calm conditions, 20°wave angle change, (b) calm conditions, 60°wave angle change, (c) stormy conditions, 20°wave angle change, (d) stormy conditions, 60°w ave angle change. Figure S4. Wave angle change (20°and 60°) sensitivity analysis for 0.030 m/m slope, (a) calm conditions, 20°wave angle change, (b) calm conditions, 60°wave angle change, (c) stormy conditions, 20°wave angle change, (d) stormy conditions, 60°w ave angle change. Figure S5. Sensitivity analyses of tidal range change from 2.56 m to 5.7 m, (a) 0.020 m/m slope, calm conditions, (b) 0.020 m/m slope, stormy conditions, (c) 0.030 m/m slope, calm conditions, (d) 0.030 m/m slope, stormy conditions. Figure S6. Nearshore volume flux under calm and stormy conditions for differing grain sizes