The morphodynamics of a swash event on an erodible beach

Abstract A high-accuracy numerical solution, coupling one-dimensional shallow water and bed-evolution equations, with, for the first time, a suspended sediment advection equation, thereby including bed and/or suspended load, is used to examine two swash events on an initially plane erodible beach: the event of Peregrine & Williams (J. Fluid Mech., vol. 440, 2001, pp. 391–399) and that of a solitary wave approaching the beach. Equations are solved by the method of characteristics, and the numerical model is verified. Full coupling of suspended load to beach change for Peregrine & Williams (J. Fluid Mech., vol. 440, 2001, pp. 391–399) yields only slightly altered swash flows, depending on beach mobility and sediment response time; a series of similar final beach change patterns results for different beach mobilities. Suspended- and bed-load transport have distinct morphodynamical signatures. For the solitary wave a backwash bore is created (Hibberd & Peregrine, J. Fluid Mech., vol. 95, 1979, pp. 323–345). This morphodynamical bore propagates offshore initially, and leads to the creation of a beach bed step (Larson & Sunamura, J. Sedimentary Petrology, vol. 63, 1993, pp. 495–500), primarily due to bed-load transport. Its height is directly related to bed-load mobility, and also depends strongly on the bed friction coefficient. The shock dynamics of this bed step is explained and illustrated. Bed- and suspended-load mobilities are quantified using field data, and an attempt is made to relate predictions to measurements of single swash events on a natural beach. Average predicted bed change magnitudes across the swash are of the order of 2 mm, with maximum bed changes of up to approximately 10 cm at the bed step.


F. Zhu and N. Dodd
Significantly, Zhu et al. (2012) observed the formation of a bed step (discontinuity in bed level), formed at the backwash bore (Hibberd & Peregrine 1979).
In this paper, we aim to bring together these strands of research by developing a mathematical model of the swash zone that, for the first time, includes shallow water hydrodynamics, fully coupled bed change, bed-load transport and suspendedload transport (concentration). The purpose is to understand the roles of both modes of sediment transport in the swash, and to see how they affect beach evolution during one swash event. In particular, we focus on the development of a beach bed step, which is a common feature of the swash (Larson & Sunamura 1993). To this end, we initially focus on the PW01 event, in order to allow comparison with earlier studies, and to examine the effects of bed and suspended load. Thereafter, we focus on the bed-step formation, by examining a swash event due to a solitary wave, which is a commonly used realistic model for a wave approaching the beach. We also use field measurements of beach change under single swash events on a natural beach to calibrate the unknown model parameters.
In § 2 we present the model equations. We then examine the PW01 event in § 3. In § 4 we present the solitary wave event, and in § 5 we estimate model parameters. Finally, we discuss our conclusions.

Governing equations
The nonlinear shallow water equations including bed shear stress described by a drag law (Soulsby 1997) are utilised to describe the flow in the swash zone, ht +ûĥx +ĥûx = 0, (2.1) wherex represents cross-shore distance (m),t is time (s),ĥ represents the water depth (m),û is a depth-averaged horizontal velocity (m s −1 ),B is the bed level (m), g is the acceleration due to gravity (m s −2 ) and c d is a dimensionless drag coefficient. In figure 1 we illustrate the situation that is considered. The bed-evolution (sediment conservation) equation including both bed and suspended load isBˆt whereq b is sediment flux due to bed load (m 2 s −1 ), which, in general, is strongly dependent onû and a weak function ofĥ,D is the dimensional deposition rate (m s −1 ) andÊ is the dimensional erosion (or entrainment) rate (m s −1 ). Here, ξ = 1/(1 − p), where p is the bed porosity. The governing equation for the transport of suspended sediment iŝ whereĉ is volumetric concentration, andĥĉ therefore represents the volume of sediment per unit area of seabed (m).
The Meyer-Peter-Müller formula (see Yalin 1977;Soulsby 1997), which is commonly used in engineering problems, is employed, whereû 0 is a representative velocity scale andÂ is a dimensional (m 2 s −1 ), empirically determined, representative bed-load sediment transport rate. Here,û crb is the threshold velocity for sediment motion as bed load (m s −1 ). It should be noted thatÂ = Aû 3 0 , where A is the equivalent dimensional constant of KD10.
We employ the entrainment model in Pritchard & Hogg (2003) and , taking the entrainment rateÊ =m e (û 2 −û 2 crs )/û 2 0 , withm e the parameter of the sediment entrainment rate (m s −1 ) of suspended load andû crs the threshold velocity for sediment motion as suspended load (m s −1 ). It should be noted that we setû crb =û crs = 0 in all the simulations in the present work for simplicity. We anticipate that this simplification will not significantly affect beach morphodynamics except, perhaps, for shingle beaches, whereon permeability effects not considered here are also significant; see appendix A.
The vector form of these four non-dimensional governing equations is (2.16) The eigenvalues of A are the roots of the polynomial equation The polynomial equation (2.17) has four roots, one of which, λ 4 ≡ u, corresponds to the transport of suspended load. The other three roots of (2.17) are denoted λ 1 , λ 2 and λ 3 , such that λ 1 λ 3 λ 2 . It should be noted that if σ = 0 and h = 0, when u > 0, we have λ 1 < λ 3 < λ 4 < λ 2 ; when u < 0, λ 1 < λ 4 < λ 3 < λ 2 . Moreover, λ 1 < 0 and λ 2 > 0 as long as σ = 0. For the solution of λ 1 , λ 2 and λ 3 we refer to Kelly & Dodd (2009. for a system with bed-load transport q = u 3 (σ = 0.01). Figure 2 illustrates the variation of λ 1,2,3 characteristics with Froude number when σ = 0.01. It should be noted that individual characteristics, λ 1,2,3 , can 'behave' as hydro-or morphodynamic characteristics, depending on the Froude number. It can be seen that there are sudden changes in these roles at critical flow conditions. From simple wave theory (Jeffrey 1976), a shock exists when the characteristics of the same family intersect, and it is a λ i shock if it is the λ i characteristics that intersect.

Numerical method
The specified time interval method of characteristics (STI MOC) (Kelly & Dodd 2009, which can resolve shocks very accurately, is used to solve (2.9), (2.10), (2.13) and (2.14) simultaneously. As the λ 1,2,3 characteristic fields associated with (2.9), (2.10) and (2.13) are genuinely nonlinear, and the λ 4 associated with (2.14) is genuinely linear, (2.9), (2.10) and (2.13) are combined to obtain total derivatives (i.e. the characteristic form) of h, u and B with respect to time, and (2.14) is used to find the total derivative of c with respect to time. Thus, (2.18) and (2.14) is written as the total derivative of c with respect to time, These four equations of (2.18) and (2.19) are solved numerically to obtain h, u, B and c in the combined load system.
2.3.1. Initial conditions Initial conditions are given for each case examined, but one general point concerns initial values for c and B. Any pre-suspended sediment must be included as an initial condition: c = c(x, 0). It should be noted that, depending on whether c(x, 0) is less than or greater than c eq , the initial B will immediately erode or accrete due to suspended load. Here, c(x, 0) = c eq in all simulations unless otherwise specified.

Seaward boundary condition
The seaward boundary is chosen so as to be far enough away from the shore that h and u at that point are uninfluenced by any wave reflected from the shore throughout the computation time. Consequently, the seaward boundary is chosen at x = −150. It should be noted that for both events examined there is therefore a region of uniform flow adjacent to this boundary, in which the flow can be specified analytically. Other dependent variables at the seaward boundary may be extrapolated from values at neighbouring points. Thus, we have h(−150, t) = h off and u(−150, t) = u off − tan α off t, where tan α off is the bed slope at x = −150, and B and c are extrapolated from neighbouring points. It should be noted that this boundary is therefore not specified based on incoming and outgoing characteristics, although characteristics at this location can be calculated straightforwardly.

Wet-dry boundary treatment
At the tip (shoreline), x = x s (t), h(x s ) = 0 and c(x s ) = c eq = u 2 s . For the solution of x s , u s = u(x s ) and B s = B(x s ) we refer the reader to Zhu & Dodd (2013).

Shock conditions
For derivations of shock conditions for mass and momentum conservation we refer to Kelly & Dodd (2010) and Zhu et al. (2012); the shock conditions are where L and R represent variables on the left and right sides of a shock, and W is the shock velocity. For the bed evolution, the change of bed level in the fixed domain [x 1 , x 2 ] is balanced by the net sediment flux into (out of) that domain, and also the net sediment settlement onto (or entrainment from) it, thus d dt Here, it is also assumed that a shock located at ζ (t) lies between x 1 and x 2 , i.e.
x 1 < ζ < x 2 . In the limit x 1 → ζ and x 2 → ζ , (2.22) becomes It is found that suspended load has no contribution to the shock condition for the bed evolution. The suspended sediment in the water column in the domain [x 1 , x 2 ] is balanced by the net suspended sediment flux inflow across the x 1 and x 2 sections and the entrainment from (settlement onto) the bed, thus d dt Equations (2.23) and (2.24) can then be written as (2.26) It should be noted that the shock condition of the bed evolution (2.25) in the combined load system is identical to that without suspended load in Kelly & Dodd (2010). For a λ 1,2,3 shock, if h L = 0, h R = 0, u L = W and u R = W, from (2.20), the shock condition for the transport of suspended load (2.26) is simplified as (2.27) This implies that the sediment concentration across a λ 1,2,3 shock is continuous, and it is determined by the genuinely linear characteristic field associated with the transport of suspended sediment. For a λ 4 shock, which is a contact wave, we have h L = h R , u L = u R , B L = B R and c L = c R .
For the shock fitting method, we refer the reader to Kelly & Dodd (2010), who have described the technique in great detail.

Model testing
The model is verified against both suspended-and bed-load-only models. In appendix B, we compare with the PH05 results, which are for suspended load, on a nearly fixed bed. Testing of the bed-load-only model is presented in Zhu et al. (2012).

PW01 swash event
In this section, we aim to elucidate the effects of bed and suspended load as represented by σ , M andẼ by modelling a swash event morphodynamically. We examine the PW01 event here. As noted earlier, although this event is not considered representative of many swash events (Guard & Baldock 2007), in terms of sediment transport it can be considered to be qualitatively similar (Pritchard 2009). Furthermore, its use allows us to verify against earlier work (see appendix B), and also to examine a swash event in which no significant interior shock formation takes place. Therefore, we consider a swash event with the same initial conditions as the PW01 swash, which then evolves morphodynamically; we refer to it hereafter as the PW01 event. It should be noted that this shoreline motion was originally derived by Shen & Meyer (1963). We refer to it here as the PW01 event because these authors, who extended the analytical solution for the shoreline motion to the whole swash and examined overtopping due to this event, provided it in a more accessible form, and made the connection with the same swash motion as the result of a dam-break problem, which interpretation we make use of as our initial condition.
The initial conditions of the PW01 swash are described by a dam-break problem over an erodible beach of initially uniform slope B = tan α off x = 0.1x (Kelly & Dodd 2010) (see figure 3). The dam is situated at x = 0 with still water on the seaward side and none on the shoreward side. The water depth behind the dam is h (x < 0, t = 0) = 1. The dam is assumed to collapse at t = 0, and the flow is dominated by gravity.  The effect ofẼ on the net sediment flux was discussed in depth by . In figure 4(e,f ), we show a variety of resulting profiles for different values ofẼ, now with a fully mobile bed. AsẼ increases there is less net movement overall. This is due to the erosion/deposition term in (2.13) (c − u 2 ) being near zero for largerẼ (consistent with PH05) for most of the swash event, because c ≈ c eq (because of the more rapid adjustment asẼ increases). Therefore, overall beach change is reduced (for fixed M) (see also figure 21 for verification against PH05 in terms of net fluxes). This is illustrated in figure 5.
Additionally, when pre-suspended sediment is present (figure 4f ) the net change increasingly favours deposition at the base of the swash. This is again consistent with PH05, and stems from the fast response time for largerẼ causing initially entrained sediment to be immediately deposited (in the lower swash).
When no pre-suspended sediment is present (figure 4e) the proportion of deposited to eroded sediment volume (deposited volume/eroded volume) for x > 0 varies from 0.06 (Ẽ = 0.03) to 0.05 (Ẽ = 0.001), peaking at approximately 0.07 (Ẽ ≈ 0.01) in between. It should be noted that  also observed a peak in the maximum positive net flux across the swash zone caused by locally entrained sediment, which essentially is a measure of the total amount of deposition in the swash zone. Moreover, the work of  (see figure 11a,c of PH05) has further indicated a peak value in the proportion of deposited to eroded sediment flux/volume. ForẼ → ∞, this ratio tends to 0, as expected.

Suspended and bed load
If we include both suspended and bed load we can now examine their effects on the beach during this event. The final bed changes in three combined load simulations are shown in figure 6, and correspond to beaches in which either bed or suspended load is dominant, or they are approximately equal. Also shown are the bed changes caused by the equivalent bed-load-only and suspended-load-only simulations, and that due only to those same components of the combined load simulation. It should be noted that, for the parameters chosen here, each mode of sediment transport has little effect on the other, and each has a distinct morphodynamical signature, and these are consistent with the results of Pritchard & Hogg (2005)    Combined change Change due to bed load Bed-load-only Change due to suspended load Suspended-load-only FIGURE 6. Bed changes in the combined load simulations and comparisons with bed-load-only and suspended-load-only simulations after one single PW01 swash: (a) bed-load dominance; (b) bed-and suspended-load approximate parity; (c) suspended-load dominance.

Swash event with shock formation
In this section, a swash event that involves shock formation in the swash zone is simulated to examine the effect of a more representative swash event on an erodible beach, and, in particular, to study the bed-step development associated with the backwash bore (Hibberd & Peregrine 1979;Zhu et al. 2012). Initially we exclude bed shear stress, and focus on the shallow water dynamics and backwash bore and bed-step developments, and consider what happens when the shoreline encounters this feature. Thereafter we introduce bed shear stress and consider the shock dynamics that contributes to the bed-step development. A swash event driven by a solitary wave is examined. The event represents a simplified but physically appropriate model of a surface gravity wave encountering a beach.

Solitary wave event
A solitary wave of height H w = 0.6 on a still water depth of h st = 1, with its crest located at x = −22 when t = 0, is considered, see figure 3. In the region x −10 the bed is flat (tan α off = 0), while for x −10 the beach is of a uniform slope tan α = 0.0667. For x −10, h(x, t = 0) = 1 − tan α(x + 10), u(x, t = 0) = 0 and B(x, t = 0) = tan α(x + 10). In the region x < −10, the water depth is determined according to the surface profile for a solitary wave given by Mei (1990) ). The water velocity is then determined by the hydrodynamic Riemann invariant along the backward characteristic: . However, the water flow is assumed not to be affected by the bed change at the initial time. Across the domain, c(

Simulation without bed shear stress
The flow structure as a solitary wave travels shorewards over an erodible beach (σ = 0.01, M = 0.001 andẼ = 0.01) simulated by the combined load model is shown in figure 7. When the solitary wave travels shorewards, it breaks and forms a shock (bore) travelling to the shore. The shock then collapses at the initial shoreline position (x = 5), and then the water climbs up the dry beach. The water velocity reaches its maximum when the shock collapses, and then decreases when the flow climbs up the dry beach. The flow (in the region x > 5) is similar to that in the PW01 swash (cf. figures 7a and 22). This is because the wave breaks at the base of the swash with little water momentum behind it, similarly to the PW01 event (see Guard & Baldock 2007).
The bore collapse at x = 5 causes considerable entrainment into the water column, which is then moved onshore. Erosion prevails at the base throughout the computation time; see figure 7(c). In the early uprush, c < c eq in the lower swash, and so sediment is transported to the upper swash (figure 7d,e). Conversely, in the later stages of the uprush, c > c eq , and suspended sediment begins to settle out and c decreases. The overall picture is similar to that depicted for moderate/smallẼ values by .
It should be noted that at the shoreline c = c eq = u 2 s and at the maximum inundation it is 0. This is because as the tip is approachedẼ/h → ∞, so that there is an instantaneous adjustment of c to c eq . For the whole backwash (for x > 5), c < c eq , and sediment is entrained into the water column and moved seawards, and the beach is eroded; see figure 7(d,e).
Bed-load transport results in erosion everywhere (not shown) (because q bx > 0 everywhere), except at the shoreline, where there is instantaneous increase (decrease) in B due to the sediment bore advancing (receding) (see Kelly  (e) instantaneous deposition/erosion distribution due to suspended load. The thick black dashed line in (e) represents the backwash bore path.
The shoreline reaches the maximum inundation x ≈ 42.6 at t ≈ 50.3. The accelerating offshore flow rapidly becomes supercritical, and, flowing into slowly retreating deeper water, results in a backwash bore (Hibberd & Peregrine 1979;Guard & Baldock 2007;Antuono & Hogg 2009; a hydraulic jump), which is also found in the morphodynamic bed-load-only simulations for the HP79 swash event (Zhu et al. 2012). The development of a backwash bore can be seen in figure 7(a,b). The strength gradually increases, with increasing differences in h, u and B on the two sides of the shock, and therefore leads to the development of a bed step coincident with the bore; this shock moves offshore and gradually slows. The shoreline eventually catches up with this shock, leading to a fully developed bed step at the base of a dry beach; see § 4.3. Across the shock u is discontinuous, and is c continuous, (2.20), (2.21) and (2.27), which results in entrainment and erosion (deposition and accretion) on the shoreward (seaward) side of the bed step as (from) suspended load (see (2.13)). This continues until the sediment concentration c on the shoreward side exceeds c eq , at which point the boundary between suspended-load erosion and accretion departs from the shock path (figure 7e).
It is the bed load that is most closely linked with the bed-step development, however, as this responds instantly to all flow changes in the vicinity of the shock. Across the shock the massively different bed-load transport rates (q b ) result in accretion immediately seaward of the shock, and development and offshore propagation of the bed step. It should be noted, however, that shoreward of the shock there is also local accretion due to convergence of bed-load transport. This is caused by corresponding initially small gradients in u (figure 7b) shoreward of the shock, which instantly affect q b (unlike q s ), leading to local deposition. Thus, the bed step both advances seawards and grows in height.
It should be noted that although suspended load can affect a moving bed discontinuity through changed erosion/accretion rates across the shock, and by modifying the beach profile generally, the bed-step height and velocity are directly linked to bed-load transport (as long as W = 0; see (2.25)).

Analysis of when the shoreline catches up with the backwash bore
In the simulation without bed shear stress, the shoreline eventually catches up with the backwash bore. As there is a sediment bore at the shoreline, this is a shock-shock collision. In this situation there are three regions, with the rightmost being dry, and in which the right shock (the shoreline) is wet-dry and the left (backwash bore/bed step) is a wet-wet shock; see figure 8.
As the two shocks come closer, the middle region gradually vanishes, and in the limiting case this region converges to one point, such that the values of dependent variables on the left and right extremities of this region become equal. Therefore, at the moment of the two shocks colliding, the middle region is assumed to disappear and the flow has only one newly formed discontinuity; however, the shock conditions 124 F. Zhu and N. Dodd are then usually not satisfied at this new discontinuity, which is therefore not stable and collapses, with the resulting states found by solution of a local Riemann problem. This is not pursued here.
As the middle region width tends to 0, h R → 0 and u R → dx s /dt = u s . As the numerical code cannot proceed when the two shocks are very close, the numerical solution is only an approximation for the case of zero middle region width. The analytical solution for the backwash bore (appendix E), of which one side is a nearly dry bed but with water of finite velocity, can be utilised to obtain the limit flow structure right before the shock collision.
From appendix E, when the shoreline approaches the backwash bore Across the sediment bore at the shoreline we have B D = B R − σ u 2 s . Thus, as the shoreline approaches the bed step Therefore, the resulting (post-collision, unstable) bed-step height will be slightly lower than the water depth on the seaward side h L . The water on the seaward side may overtop the bed step, depending on the velocity and also bed mobility.

The inclusion of bed shear stress
The contour plots for the solitary wave simulation with bed shear stress (σ = 0.01, M = 0.001,Ẽ = 0.01 and c d = 0.01) are shown in figure 9. It should be noted that the shoreline does not recede in the backwash (Antuono, Soldini & Brocchini 2012;Zhu & Dodd 2013), because of bed friction, and the swash zone is always wet once wetted. It should be noted that the corresponding swash period is much longer than for c d = 0. The flow structure in the uprush is similar to that in the simulation without bed shear stress, with an overall deeper flow and smaller velocity. However, bed shear stress greatly reduces the seaward velocity of the backwash flow, and the backwash duration exceeds that of the uprush. The bed shear stress also modifies the flow structure in the backwash, and more shocks are formed in the backwash. The shock paths (including inception and termination points) are illustrated in figure 9(a-c). It should be noted that the first shock to form (a λ 1 shock) does so at about t = 23.5. It quickly propagates offshore, and has no influence on the subsequent bed-step development, so we ignore it hereafter. It should be noted that the values of c in figure 9(d) are considerably reduced compared with those in figure 7(d), especially during the backwash, consistent with reduced backwash velocities.
In the backwash, a bore (λ 3 shock) develops (the backwash bore), which initially travels seawards (see figure 9). This is a robust feature of the backwash (see also figure 7a). At t ≈ 61.4, another weaker λ 3 shock collides with the backwash bore. Here, we treat the shock collision by neglecting the weaker shock when it is in close proximity to the backwash bore. This approach shows good agreement with the idealised Riemann solution. The collision increases the backwash bore (λ 3 ) shock strength slightly; see figure 10. Thereafter, the backwash bore slows, and decreases in strength (see figure 10).
At t ≈ 67.7, a λ 2 shock develops on the shoreward side of the backwash bore/bed step, the strength of which rapidly increases (see figure 10), and which propagates shorewards, as the flow in the later backwash diminishes. Thereafter it diminishes and (d) c. In (a-c) we show the shock paths of the λ 3 and λ 2 shocks, with their points of inception and termination. Dotted line: λ 1 shock; dot-dashed line: first λ 3 (weak) shock; dashed line: second λ 3 shock (backwash bore); solid line: λ 2 shock. slows, in effect re-establishing a 'shoreline' boundary between the sea and water still draining back. It should be noted that, as with all morphodynamical shocks, this λ 2 shock is accompanied by a bed discontinuity, which travels with the shock, but which is much smaller than the bed step associated with the original λ 3 backwash bore/bed step; see figure 10.
At t ≈ 73.9, the backwash bore/bed step (λ 3 shock) is brought to rest: W = 0. This is the point of reversal. From the shock conditions, a bed discontinuity with W = 0 exists only when u L = u R = 0. Furthermore, h L + B L = h R + B R when W = 0, indicating the continuous surface level across the bed step (λ 3 shock) at the point of stationary state. Thus, there is at this instant only the bed step, and an associated change in depth, with other quantities being continuous. Thereafter u L and u R increase, but the backwash bore (λ 3 shock) vanishes, because the characteristics on either side of the bed step diverge, thus allowing the shock conditions only to be satisfied by a nonphysical shock (i.e., one with diverging characteristics). Mathematically, the structure at the former bed-step position after reversal could be derived from an idealised dambreak problem, and (seaward to shoreward) is a λ 1 shock, a λ 3 rarefaction fan and a λ 2 shock. The rarefaction fan is almost stationary, and the change in variables across it is much larger than that across the two (weak) shocks (which are not shown here). Thus, A E represents the inception of a shock.
the previous bed discontinuity acquires a continuous structure; however, the observed bed-step profile is little changed. This process of the original seaward moving λ 3 shock (backwash bore/bed step) slowing, and the creation of a fast-moving mainly hydrodynamical λ 2 shock, is equivalent to the reversal of a hydraulic jump on an immobile bed but here results in a bed step (Larson & Sunamura 1993) feature being left, almost inert, on the lower beachface as the hydraulic jump reverses. Finally, it should be noted that all three shocks are formed at approximately critical conditions (F r ≈ −1) (figure 10d). It is for these conditions that the dispersion curves for the λ 2 and λ 3 waves (figure 2) change slope rapidly, implying that relatively small changes in flow may result in convergence of characteristics. Further, it should be noted that as an initially small shock passes through critical conditions it is likely to grow, for the same reason (it should be noted that the definition of shock strength here is the jump in characteristic slope across the shock); see also figure 10(a,d). It should be noted that the weak λ 3 shock does not pass through critical conditions. Finally, in figure 11 we show the characteristic fields associated with the shock dynamics described above. The characteristic convergence and shock formation can clearly be seen.

4.5.
Final beach change and bed-step development The beach changes due to different σ , and M when the shoreline catches up with the backwash bore in the simulations without bed shear stress are shown in figure 12. For all the cases examined, a bed step forms, the height of which decreases as σ decreases. The bed-step height decreases slightly as M increases, due almost entirely to the increased bed level on the seaward side of the bed step, which itself is due to greater sediment deposition for increased M. Thus, bed load (σ ) dictates crest level (figure 12), but suspended load (M) can affect height. This is also shown by figure 13. The effect of M on the step height is not large compared with that of σ , but becomes more pronounced for small σ . It should be noted, however, that these heights refer to the size of the jump in B (B R − B L ); the overall prominence of the step is dictated by σ ( figure 13). One should note also the increased erosion shoreward of the step for larger M (figure 12).  To examine the importance of the swash event to the bed-step generation, the final water surface profiles and bed changes for both the present solitary wave simulation and the HP79 swash (see Zhu et al. 2012) (both for c d = 0) are illustrated in figure 14. It should be noted that both events have the same wave height, but that the HP79 event is essentially an adjustment in water level, so that the backwash (and therefore backwash velocities) is (are) much reduced compared with the uprush. It should be noted also that the bed steps for the solitary wave case including friction (figure 15), in which backwash velocities are significantly reduced and no longer dependent solely on maximum run-up, are consistent with those for the frictionless HP79 case. The bed-step heights in both swashes are close to but slightly smaller than the water depths due to the sediment bore at the tip, consistent with the analysis in § 4.3.
In figure 15, we show final bed changes for various c d values. (Profiles are deemed 'final' when h R 0.006, where here h R is the water depth on the shoreward side of the shock nearest to the shoreline.) Increasing bed friction strongly influences the bed-step development, both the step size and its location, by reducing backwash velocities, and thereby the λ 3 shock strength and its accompanying bed-step size. Finally, it should be noted that the bed step shown in figure 15 (except that for c d = 0) is no longer a discontinuity, but in fact a rapid variation in B.

Estimation of M and σ from field experiments
Finally, we note that it is difficult to evaluate M for real beaches (Ẽ may more straightforwardly be evaluated from settling velocities). In an attempt to do this we present figure 16, in which we show a contour plot of net onshore flux of suspended sediment entrained in the uprush only at a location in the mid-swash (x = 13) of the solitary wave swash event (Q en,up (13); see appendix C for the definition of net flux Q) with bed shear stress (c d = 0.01) as a function of M andẼ. This position is roughly equivalent to that of the sediment traps in Masselink & Hughes (1998) and Hughes, Masselink & Brander (1997). Both studies, which were for grain sizes and beach slopes of, respectively, 0.5 mm and 0.14, and 0.3 mm and 0.12, yielded representative average onshore fluxes in moderate wave conditions of ∼30 kg m −1 over one swash uprush. Thus, a grain size of 0.4 mm (w s ≈ 0.05 m s −1 ⇒Ẽ ≈ 0.02 ifĥ 0 = 1 m) corresponds to a value of M ≈ 0.001.
As with determining M, determining σ is an uncertain process, but its effect can be quantified similarly to that of M in the variation of the net uprush bed-load flux in the mid-swash (where t in (t de ) is the time of inundation (denudation)); see figure 17. A net flux of 30 kg m −1 due to bed load thus corresponds to σ ≈ 0.01. Further, it should be noted that figure 6 implies that assuming that suspended-and bed-load transport do not significantly affect each other, as done in figures 16 and 17, is reasonable. Figures 16 and 17 thus allow a determination of M and σ based on a representative bore-driven swash event on a sandy beach. Although we can obtain representative midswash uprush fluxes from the field measurements of Masselink & Hughes (1998) and Hughes et al. (1997), we do not know the proportion of bed and suspended load in these experiments, and the above estimates implicitly assume that all transport is by one mode only. The consensus appears to be that both modes are important in the uprush (Horn & Mason 1994;Masselink & Hughes 1998), in which case we may  Bed mobility parameters M and σ , and initial concentrations for four scenarios in which the net mid-swash onshore flux is 30 kg m −1 . Also shown are the height of the resulting bed step the volume of sand (V) in the step, the average bed change magnitude from the initial shoreline position (x = 5) to the maximum inundation (x max ) and that from x = −5 tox max (therefore, that including the bed step), and the net volumetric sediment transport relative to the initial shoreline position.

Proportions of bed and suspended load
We use figures 16 and 17 to allocate a proportion of the nominal 30 kg m −1 to each mode. These allocations are summarised in table 1, and correspond roughly to a scenario in which bed and suspended load are both significant (test 1), which therefore comprises our best estimate of reality, and two further ones in which bed (test 2) or suspended (test 3) load dominates. Finally, in test 4, we consider a case in which pre-suspended load dominates. Therefore, for test 4, we reduce the local bed-and suspended-loadQ values (and therefore M and σ ) and then impose c at the time of bore collapse at the shore: t = t c . This is to reproduce the large local concentrations that might be expected due to bore turbulence, which is not present in the mathematical model. The approach here is to impose c(x, t = 0) = c eq (x, 0) and c(x, t c ) = nc eq (x, t c ), and to determine n such that the net total mid-swash onshore flux is 30 kg m −1 . This results in n = 1.8. It should be noted that pre-suspended load exists in tests 1-3 too, but only at equilibrium values (at t = 0). Therefore, much of the sediment initially present falls out of suspension prior to bore collapse; see § 4.1. The volume of sand in each bed step (V) is also presented in table 1, and is calculated from the base on the left side B = B L to the position on the right side such that B R = B L , see figure 19(b).

Concentrations and final bed profiles
The contour plots ofĉ for four test cases are shown in figure 18. The suspendedload dominant case test 3 has the largest concentration with a maximum of ∼0.01. In the uprush the effect of the pre-suspended sediment can still be seen, but not in the backwash ( figure 18d). The corresponding final beach profiles for these four tests are shown in figure 19. In figure 19(a) we show bed changes at the point of λ 3 shock reversal (W = 0), at which point the bed step achieves peak magnitude. At this instant the bed-step height is directly proportional to σ , as noted earlier, as is the volume (see table 1). Thereafter, the degree to which the bed step is modified depends on M primarily, with the largest reworkings in test 3 (M = 0.0009) and test 1 (M = 0.00066); see figure 19(a).
Ifĥ 0 = 1 m, we obtain average bed change magnitudes in the region shoreward of the initial shoreline position of about 2 mm, with moderately larger values including the bed step. Masselink & Hughes (1998) and Hughes et al. (1997) do not record the bed change due to individual events, but Blenkinsopp et al. (2011), on a beach with a grain diameter of 0.4 mm, but of slope 0.067, record most (∼60, 75, 95 % moving from seaward to shoreward in the swash) individual swash events with a zero (or at least non-measurable) net beach level change at three cross-shore locations across the swash zone. These negligible changes would appear to include beach changes similar to the average magnitudes calculated here, from the initial shoreline position (figure 7 of Blenkinsopp et al. (2011)). Nonetheless, there is a significant proportion that has non-zero net change (both positive and negative) up to an occasional extreme of ∼4 cm. Based on this comparison the changes recorded here are a little larger than those observed, particularly at the beach step, but only moderately so. It is not clear where the measurements of bed-level change were made relative to where a bed step might form, but it seems likely that some of them might have been made at such a location, because the measurements of Blenkinsopp et al. (2011) were made over tidal cycles.
Finally, it should be noted that the values ofĉ in figure 18 are consistent with values measured in some field experiments .

Concluding remarks
A mathematical model in which, for the first time, the 1D shallow water equations (with bed shear stress) are fully coupled both to an Exner equation and a concentration advection equation is presented.
Numerical simulations of one single PW01 swash reveal that the sediment entrainment rate parameter M controls the amount of erosion/deposition caused by suspended load, and that the resulting beach change patterns are similar, with amplitude ∝ M (ifẼ, the parameter governing the sediment response rate, is constant). Simulations also reveal that coupling with suspended load has only a minor feedback on the flow itself. Furthermore, bed-and suspended-load transport do not significantly affect each other.
Simulations of a single solitary wave swash event reveal some important features of swash flow. A bed step associated with a backwash bore (a λ 3 shock in the present system) is formed in all simulations, with height proportional to σ (bed-load mobility) and inversely proportional to c d (Chezy coefficient). The results show that suspended load can also affect bed-step height, although only slightly, by changing the bed level on the seaward side due to sediment deposition. This is consistent with the analytical result that suspended load has no direct effect on the shock conditions. However, subsequent sediment entrainment after flow reversal, when the bed step is transformed into a (steep) rarefaction fan, can modify its height.
The bed step grows as the backwash bore gradually slows down, and achieves a maximum amplitude at a stationary state with u L = u R = W = 0. After this point the backwash bore (i.e. shock) vanishes and the flow starts to move shorewards.
The previous bed discontinuity acquires a continuous structure, although the observed bed-step profile is little changed.
Near to where the point of reversal of the bed step is encountered, but before it, a shoreward propagating λ 2 shock forms (i.e. near the end of the swash event). It grows rapidly and is the main mechanism for re-establishing the shoreline as the swash motions decay. Both λ 3 and λ 2 shocks form near to and grow rapidly as they pass through critical conditions. This process of the original seaward moving λ 3 shock (backwash bore/bed step) slowing, and the creation of a fast-moving mainly hydrodynamical λ 2 shock, is equivalent to the reversal of a hydraulic jump on an immobile bed, but here results in a bed-step feature being left on the lower beachface.
There are a number of limitations to the present study. The most obvious is the swash event itself and the fact that M and σ as evaluated here are therefore empirical values. It is likely that different swash events will yield rather different pictures of erosion and accretion in the region. Furthermore, it is noted that the swash events that move the largest amounts of sediment, at least in some studies and on some beaches, are usually those that include one or more interactions (Blenkinsopp et al. 2011). Therefore, some circumspection is required in interpreting the present findings. Nonetheless, a solitary wave is a robust and widely accepted model for a wave approaching the shoreline for steeper beach slopes. Furthermore, the use of uprush sediment transport only as a characterisation of bed mobility is more robust than that of net transport. The value chosen (30 kg m −1 ) is consistent with a number of field studies (Hughes et al. 1997;Masselink & Hughes 1998;Blenkinsopp et al. 2011), and such an event might be characterised as a moderately large but not exceptional swash event. It should be noted also that in reality M, σ ,Ẽ and c d are all related to grain size.
The present model also neglects bore turbulence, which, as noted earlier, would entrain and suspend more sediment if included. It is therefore possible that M might be overestimated here to compensate for this. However, beach mobility will also affect entrainment by bore turbulence, so its basic effect is robust.
As mentioned earlier, the effect of a threshold of motion for sediment is considered in appendix A. It is not considered to be significant.
The formation of a beach (bed) step is one of the most interesting features of the present study. This is a realistic beach feature (Larson & Sunamura 1993;Masselink et al. 2010). In the field it has been reported to reach about 0.5 m in height (Masselink et al. 2010), albeit on a more permeable beach with grain size ∼5 mm. The nearly vertical slope predicted here is, of course, unrealistic, but no significance should be read into this because downslope diffusion (avalanching) under gravity is excluded, to allow an understanding of the shock dynamics from which the step forms, and becomes (relatively) inert. In the field angles of ∼20 • are typical (Larson & Sunamura 1993), although this slope is presumably proportional to grain size, with coarser grained beaches on which steeper slopes will form also being more subject to permeability effects.
It is suggested that future studies should investigate the role of different types of, and multiple, swash events, not least because the backwash bore and bed step are likely also to be important in entrainment . Masselink et al. (2010) also report that the berm height (the berm being the depositional region at the top of the swash zone) could be controlled by the step size. In order to investigate the effect of suspended load on the swash flow, we show the comparison of contours of h and u in the suspended-load-only simulations (σ = 0, M = 0.001/0.005 andẼ = 0.001) with the equivalent fixed-bed case (σ = 0 and M = 0), i.e. the PW01 solution, in figure 22. The flow structure of the suspended-load-only simulation with M = 0.001 is very close to the PW01 solution, although the final bed profile is changed to a certain extent (see figure 4a); the maximum net beach change ≈0.023. For M = 0.005, the flow is changed to a greater extent due to the larger bed change (see figure 4a) (the maximum net beach change is ≈0.11). However, the maximum inundation is changed little from that in the PW01 solution. Even though the final bed change for M = 0.005 is comparable with that for 0.01 < σ < 0.0654, the flow for the suspended-load-only simulation is little changed in comparison with the bed-load-only simulation. The smaller effect of suspended load on the swash flow indicates the lesser importance of fully coupling for the suspended-load-only simulation.
Appendix E. Shock relation when one side of the shock is a (nearly) dry bed for the combined load system One special case occurs when one side of the shock is a (nearly) dry bed. Here, we examine the shock behaviour when its right side is becoming dry with h R → 0. The case of h L > h R is considered.
From (2.20), we have where m represents water mass flux across a shock front. As h L > h R , we have |u R − W| > |u L − W|. If m > 0, we have W < u L < u R , and if m < 0, u R < u L < W. From (2.21), At the limit of h R → 0, m → 0 and the right-hand side of (E 2) approaches 0. Thus, it is possible that h R + B R − h L − B L → 0 or (and) h L + h R → 0. However, which term approaches 0 at the limit depends on the sign of m and W. It is therefore classified into the following four cases according to m and W to find the solution to the shock adjacent to a nearly dry bed.
Hence, with h L > 0, we have h L + h R > 0, and when h R → 0, i.e. m → 0, it is (h R + B R − h L − B L ) → 0 such that (E 2) is satisfied. Thus, at the limit h R → 0 we have h L → h R + B R − B L > 0. (ii) m > 0 and W < 0 When m > 0 and W < 0, we also have u R − u L > 0, but B R − B L = σ ((u 3 R − u 3 L )/W) < 0. Similarly, we still have When h R → 0, h R + B R − B L < 0, and as h L 0, h R + B R − h L − B L < 0. Thus, it has to be h L + h R → 0 such that (E 2) can be satisfied. This gives h L → 0 when h R → 0.
In summary, the shock solution at the limit can be obtained according to the travelling direction of the shock and that of the water mass across the shock front. In cases (i) and (iv), h L → h R + B R − B L → B R − B L when h R → 0; meanwhile, in cases (ii) and (iii), h L → 0 as h R → 0.
From (E 1), W = u L + ((h R (W − u R ))/h L ). As h R < h L , h R /h L → 0 when h R → 0, and we have W → u L regardless of whether h L → 0 or not. The two shock conditions (2.20) and (2.21) have been simplified into two relations when h R → 0. Furthermore, (2.25) is also used to solve the shock. Thus, the system is determined. When the signs of m and W are determined, the shock solution is unique for a shock with nearly dry bed on its right side. According to the criterion for a physical shock that the characteristics must converge, cases (i) and (ii) are not physical.
For the backwash bore, m < 0, W < 0 and u R − u L < 0, and it corresponds to case (iv). Thus, we have h L → B R − B L > 0 as h R → 0.
Finally, we also note that using (2.20), (2.21) and (2.25) we may write where an overbar denotes a simple average (e.g. h = (h R + h L )/2), which gives us a direct relationship between h R − h L and B R − B L .