A PROCESS-BASED NUMERICAL MODEL OF SHOREFACE PROFILE EVOLUTION

The paper presents an effective integrated framework and process-based dynamics in the form of new software that successfully models the whole shore-face zone seamlessly from relative deep water to and across the surf zone to the beach. As an advance on other available models, it simulates both erosion with bar development and beach recovery after erosion events, with shoreward bar migration and beach accretion. It incorporates analytical descriptors of the hydrodynamic and morphological processes, including simplified empirical relationships that describe the forcing of waves and currents and the sedimentation responses that change the profile. The model represents the surf zone wave roller and its effects in driving undertow. The simple bed load formulation of Ribberink and Al-Salem (1990) based on 𝑢 3 ��� is used for sand transport, incorporating the effects of wave asymmetry, acceleration skewness, boundary layer streaming and Stokes drift, with the coefficient increased to account for suspension transport in the surf zone. Importantly, bar migration rather than diffusion is achieved by management of acceleration skewness as the wave passes over the bar crest. It is shown by validation to simulate profile evolution satisfactorily in a manner and at time-scales consistent with measured data, for both small scale laboratory and long term prototype modal and storm conditions. The model caters for water level variations due to tide, storm surge and/or sea level rise. It may be used to assess beach nourishment deposition at the beach berm, lower shore-face or across the surf zone.


INTRODUCTION
It remains a challenge to model numerically profile evolution across the shore-face from deep water to the shoreline with stability and accuracy.At long time scales, attempts to date typically utilize equilibrium profile shapes or bed slopes, with profile change based on the degree of disequilibrium and wave/current forcing, primarily for the surf zone in erosion events.Simulation of bar dynamics, particularly shoreward bar migration and beach recovery after erosion events, is generally poor due to uncontrolled morphology feedbacks.
Here we present an effective integrated framework and process-based dynamics in the form of a new software package that successfully models the whole shore-face zone seamlessly from relative deep water to the beach.It incorporates analytical procedures, including some simplified empirical relationships that describe the forcing of waves and currents and the sedimentation responses that change the profile.It is shown to simulate profile evolution satisfactorily, with shoreward as well as seaward sand transport, bar development and migration and beach erosion and accretion, in a manner and at time-scales consistent with measured data, for both small scale laboratory and long term prototype modal and storm conditions.The model caters for water level variations due to tide, storm surge and/or sea level rise.It may be used to assess beach nourishment deposition at the beach berm, lower shore-face or across the surf zone.
The model aim is to simulate time-dependent evolution of 2DH beach profiles under constant or time-varying wave and water level conditions over unlimited durations.The model framework contains many process components that function interactively to simulate the time-dependent wave conditions, hydrodynamics and sedimentation behaviour at each time-step.The model provides for input of: • The initial profile, defined either as a laboratory plane slope or by measured x-z data for prototype cases; • Sediment d 50 (m) and water temperature, which are used to estimate the settling velocity w (m/s); • The swash slope, berm level and dune height • Waves as H, T time series, (H s and T p for irregular waves) and water depth of wave data; • Water level time series for tides, surges and/or SLR.
The model utilises the 1-D Exner relationship for bed level change across the horizontal x domain from the seaward shore-face boundary to the shoreline.It may be applied over time-scales of hours to months/years at time-steps of (order) 0.1 to 1 second for laboratory and prototype scales.

Wave Propagation and Hydrodynamics
Waves are propafated initially using linear theory and converted empirically to finite amplitude heights to determine wave breaking conditions and to force hydrodynamic and sedimentation responses.Wave breaking is determined on the basis of breaker index () = (    ⁄ , ()), which varies with bed slope (s) such that it approaches about 0.5 for s=0.0.

Surf Zone Roller and Undertow
Previous research has identified the role of the roller in forcing undertow (Svendsen, 1984).The roller is generated by wave breaking and dissipates as its energy is lost to turbulence and other factors.Thus, surf zone waves are considered in two parts: an underlying organized breaking or unbroken wave; and a surface roller (Figure 1).Attempts have been made previously to describe the form and size of the roller (Svendsen, 1984;Dally and Brown, 1995).Here it is assumed to have an ad-hoc empirical shape approximating its generation proportional to the breaking of the underlying wave and progressive decay that occurs both during and after cessation of breaking.The roller profile adopted in the model (Equation 1) is conceptually consistent with observed processes and has been derived to yield the most suitable results when used to drive undertow.It remains a matter for further research to incorporate a more physicsbased descriptor.
Its shape incorporates a distortion factor (D) based on that of Vellinga (1983) that provides for the effect of wave height scale and sediment settling velocity (w) on profile distortion within the surf zone.Here, D is applied to the forcing hydrodynamics to achieve profile distortion.The Vellinga formulation for D is modified to include the effect of wave period via steepness (    ⁄ ), as in Equation ( 2 (2) It is assumed that the radiation stress gradient in the underlying organized wave forces a setup profile across the surf zone while the undertow is determined as the friction velocity associated with the potential setup gradient forced by the total roller profile gradient, in accordance with Equation 3.. (3)

Sediment Transport
Sediment transport q(x,t) is the sum of shoreward wave boundary layer sheet flow transport due to asymmetry, streaming and acceleration skewness under the unbroken and underlying surf zone waves; and seaward transport driven by undertow and Stokes drift.The balance of these determines its net rate and direction, erosion or accretion, bar development and profile evolution.We adapt the simple Bailard (1981) relationship for wave boundary layer bed load  =  1  3 ��� (solid m 3 /m/s), with the Ribberink & Al Salem (1990) value  1 = 0.00018 (s 2 /m) for total load on a horizontal bed.Here,  3 ��� due to wave asymmetry is modified to include acceleration skewness, boundary layer streaming, undertow and Stokes drift.The coefficient c 1 is increased across the surf zone to account for suspended load under turbulent wave breaking and undertow.
Local bed slope (s) moderates the shoreward horizontal bed load transport in accordance with (, ) = (, 0)�1 −    () ⁄ � (Patterson, 2013;Patterson & Nielsen, 2016), in which S eq (x) is the local wave and water depth dependent equilibrium bed slope determined as () 1 , where  = and D 1 is a distortion factor = �  , , � relative to a reference condition.Following the methodology of Vellinga (1983), we establish a reference equilibrium profile based on the measured Gold Coast cyclone erosion bar and shore-face shape extending seaward of the depth of closure to 20m depth, where considerable transport occurs but the profile shape is maintained (Figure 2).Adopting reference cyclone waves H s =6m and T p =11s yields (  ) = 0.0082  0.31 .Analysis of a range of laboratory and prototype data for other wave and sediment scales yields: For regular waves, it is assumed that H and T are equivalent to H s and T p .
Effective Near Bed   ����  Sand transport utilizes an effective near bed  3 ���  that comprises contributions from wave asymmetry and near-bed steady currents: Stokes drift; surf zone undertow; and boundary layer streaming.Wave asymmetry is considered in terms of both steady finite amplitude waves of high order Fourier shape, which are asymmetric about the crest with velocity skewness (S u )>0 and acceleration skewness (S a )=0, and waves propagating on a sloping bed, which exhibit S a >0 (Figure 3).Near bed  3 ���  associated with Fourier asymmetry is determined as: For waves propagating on a sloping bed, Boussinesq tests provide  �() that includes non-zero orbital current acceleration, from which we derive  3  ��� that is expressed in the form of Equation ( 6) after Nielsen and Callaghan (2003) and Nielsen (2009), with ≈45 o .
For computational efficiency, the values of  3 ���  and  3 ���  at the bed are estimated using empirical relationships derived from a wide range of steady high order Fourier analyses (Fenton, 1988) and shoaling Boussinesq model (FUNWAVE) cases.Relationships with and without acceleration skewness have been derived, normalized with respect to the finite amplitude  2 ���  1.5 , as functions of the Ursell Number ( =  2 /ℎ 3 ) (Figure 4).For simplicity,  2 ���  is calculated from the linear  2 ���  as shown in Figure 5 and  is calculated using the finite amplitude H fin , derived from the linear theory H lin using the simple analytical relationship given by Nielsen (2009) based on data from Sato et al (1992): Importantly, we find that the mean acceleration �   � � 3 ������������ may vary with bed slope and in wave breaking.Boussinesq tests show that it diminishes in non-breaking waves passing over bar crests, leading to only Fourier asymmetry on the lee side.This affects the spatial distribution of q x across the bar crest.We incorporate an ad-hoc transition from  3 ���  to  3 ���  during that process, which significantly improves simulation of shoreward bar migration rather than diffusion.We assume also that acceleration skewness reduces in proportion to the extent of wave breaking across the surf zone.In the presence of an additional steady near-bed current  �,  3 ��� is approximated as: Parameters t c , t t , U c and U t are crest and trough durations and peak currents respectively.The expression in square brackets of Equation ( 8) represents the incremental contribution of the near-bed current to the total  3 ��� and is used with coefficient c 1 to calculate the incremental contributions to q x of boundary layer streaming, Stokes drift and undertow.The exception to this is that c 1 is factored up in the turbulent surf zone to provide approximately for additional suspension transport.

Ripple Regime Transport
Cross-shore sand transport is affected by bed ripples and may be directed seaward, opposite to the forcing asymmetry, in the presence of strong vortex ripples.Nielsen (1992) showed that ripples form at mobility number () values less than about 300, while sheet flow occurs at higher values.

MODEL VALIDATION Laboratory Wave Flume Cases
The model has been assessed by correlation with various laboratory profile evolution experiments as well as field data conditions, as an assessment of overall performance and, particularly, the distortion factors D and D 1 in equations ( 2) and ( 4) respectively.Figure 7 shows profiles and q x modelled over 5 hours correlated with the measured 5 hour profile for Experiment SE1 of Atkinson (2018).The initial

Large Wave Tank Regular Waves
A series of profile evolution tests were undertaken in the Large Wave Tank at the CERC Waterways Experiment Station in March 1956-57, as collated in Kraus & Larson (1988).

Measured
Figure 9 shows correlation of measured and modelled q x and profiles at 20 hour intervals through the 100 hour experiment.While the behaviour in the flume is not comprehensively reproduced, the model simulates the broad spatial and temporal q x and profile evolution reasonably well in both magnitude and shape, including the extent and nature of erosion of the beach and bar development.

LONG TERM SIMULATIONS OF FIELD CONDITIONS Erosion and Recovery
The capability of the model for both erosion by storm waves, with development of multiple bars, and beach recovery with shoreward sand transport and migration of the residual storm bars is illustrated in Figure 10.In the case shown, an initial profile based on typical modal conditions at Gold Coast, Australia is subjected to a hypothetical one year wave time series beginning with a major cyclone event (H s =6m) and with subsequent modal conditions (H s =0.5m-3.0m)as illustrated in the top panel.The modelled erosion profile immediately following the cyclone event is shown in the second panel (b), with development of three bars.Subsequent evolution with shoreward bar migration and recovery of the beach berm is shown in the third panel (c).The bottom panel (d) shows the change of beach width (  −   ) over the 1 year simulation, with initial erosion followed by gradual recovery.Lower Shore-face Sand Transport Surveys since 1966 of northern Gold Coast profile transects have identified a disequilibrium lobe at the lower shore-face caused by sand deposition associated with migration of the Nerang River mouth, where progressive long term evolution towards equilibrium has been used to calculate long term average rates of shoreward sand transport (Patterson, 2013;Patterson and Nielsen, 2016).Recent surveys in 2007, 2010 and 2012, together with availability of wave data for that 5 year period, provide a basis for assessing the capability of the model to determine the average q x (x) across the lower shoreface.Transect 73 (Figure 11) has been modelled using the wave time series shown in Figure 11   Comparison of the modelled and measured rates is shown in Figure 12, indicating very good agreement, with progressive increase in q x associated with shallowing depth and relatively low bed slope.The rapid decrease in q x around 1900m<x<2100m marks the development and subsequent shoreward migration of a major storm bar during 2009, which was not surveyed but was shown to have occurred by the model.the shore-face.net effect of the nourishment is determined relative to the beach width variation without nourishment.
As expected, both cases trend towards a net increase in beach width of about 10 metres.The response to berm placement is rapid initially while the shoreward sand transport and net gain at the beach for shore-face placement is gradual.The model thus provides a useful tool for assessing costbenefit of nourishment placement alternatives.

DISCUSSION
The model presented is effective in simulating erosion and bar development during storm events as well as bar dynamics including shoreward migration and beach recovery, at time-scales consistent with measured data, for both small scale laboratory and long term prototype conditions.. It successfully models the whole shore-face zone seamlessly from relative deep water to the beach.Its effectiveness derives from the model framework that incorporates simplified analytical procedures, including some empirical relationships, that comprehensively describe the forcing of waves and currents and the sedimentation responses.
These include provision for forcing of surf zone undertow by the wave roller, which is described here by an ad hoc relationship that mimics its growth during wave breaking and subsequent decay, thereby achieving the known delay in setup and undertow from the breakpoint.While it yields good results over the range of conditions assessed, it is unclear how broadly applicable it is.Ongoing  research is being undertaken to incorporate a more physics-based definition of the roller and undertow based on Dally and Brown (1995).
The model as presented extends the bed load formulation of Ribberink and Al-Salem (1990) for sand transport   =  1  3 ��� , with c 1 =0.00018 for m 3 /m/s.The near-bed  3 ���  includes provisions for both wave asymmetry and uniform currents due to boundary layer streaming, undertow and Stokes drift, The coefficient c 1 is factored up in the surf zone to cater for turbulence and suspension there.A more comprehensive approach to suspended load transport is the subject of further research.Factors are also included to cater for the effects on bed load of bed slope and rippled bed forms.
Importantly, we find that variations in acceleration skewness affecting  3 ���  with wave propagation over bars and in wave breaking play an important role in successfully simulating bar migration.Boussinesq tests show that it diminishes in non-breaking waves passing over bar crests, leading to only velocity asymmetry on the lee side.We incorporate an ad-hoc approximation to that process, which significantly improves simulation of shoreward bar migration rather than diffusion.
The model applies to long term prototype simulations involving varying waves and water levels, for which the parameters and coefficients set within the code have been shown to apply satisfactorily, with good results that correlate well with measured profile behavior.It has been shown to simulate beach and profile evolution and may be applied to assessments of sea level rise impact and beach nourishment responses.

Figure 2 :
Figure 2: Adopted Gold Coast equilibrium cyclone profile

Figure 3 :
Figure 3: Typical non-linear Fourier and Boussinesq orbital current patterns
suggested that the bed load relationship of Ribberink and Al-Salem could be moderated with a factor () and utilized rippled bed laboratory data ofRoelvink (1988) to develop the relationship () = 1 − [((500 − )) ⁄ 450] 4.5 shown in Figure6to account for reduction in q x .

Figure 6 :
Figure 6: Factor () applied to account for rippled bed effect in reducing qx
Figure 8 (top) shows the measured evolution over 100 hours of an initially planar bed of slope 1:15 subject to regular waves of H=1.52m and T=3.75s (Case 500).The sand size was d 50 =0.22mm.This represents an unusually steep wave (     � = 0.075) used to assess the model performance.Figure 8 (bottom) shows the progressive modelled profiles at 10 hour intervals together with comparison of the modelled and measured profiles at 100 hours.

Figure 8 :
Figure 8: Modelled large wave tank experiment Case 500 reported in Kraus & Larson (1988).Measured profiles over 100 hours (top) and modelled evolution with measured profile at 100 hours (bottom)

Figure 9 :
Figure 9: Modelled large wave tank experiment Case 500 reported in Kraus & Larson (1988).Measured versus modelled qx (left) and profiles (right) each 20 hours over the 100 hour test.

Figure 10 :Hs
Figure 10: Modelled prototype erosion and accretion: (a) 1 year wave Hs time series; initial cyclone peak to 6m; (b) storm bar profile; (c) progressive recovery profiles; (d) beach width erosion and recovery over 1 year.

Figure 11 :
Figure 11: Location of the wave recorder and Transect 73 at northern Gold Coast (left) and 5 year wave height time series (right).Note the occurrence of major storms during the first half of 2009..

Figure 14 :
Figure 14: Modelled profile responses to nourishment to the beach berm (top) and shore-face (centre), with effect on beach width (movement of +0.5m contour) over 1 year (bottom).
planar slope was 1:10.The waves were irregular with H s =0.15m and T p =1.3s.The sand had a reported size d 50 ≈0.3mmwith fall velocity w≈0.03m/s.The close agreement indicates good spatial and temporal model simulation of both q x and profile change.