Surface-groundwater ﬂow numerical model for barrier beach with exﬁltration incorporated bottom boundary layer model

A surface-groundwater ﬂow model is developed for the swash ﬂow on a barrier beach. The non-linear shallow water equations are used to simulate the surface ﬂow. Laplace’s equation is used to describe the groundwater ﬂow and is solved using the Boundary Integral Equation Method to provide potential heads and normal potential derivatives at and within the boundaries of the barrier. An exﬁltration incorporated bottom boundary layer sub-model is used to obtain bed shear stress. The groundwater model is veriﬁed against the numerical test results in Kazemzadeh-Parsi and Daneshmand (2012) for the groundwater ﬂow through a rectangular dam. The coupled surface-groundwater model is validated against the prototype-scale BARDEX II experimental results (Turner et al., 2016). The steady-state groundwater ﬂow comparisons show excellent agreement in phreatic surfaces. The comparisons of groundwater ﬂow under the action of waves show reasonably good agreement with experimental results in phreatic surfaces. The simulated time averaged pore velocities for the runs with and without waves are in satisfactory agreement with experimental results in general, and certain discrepancies are observed near the beach side. The bed shear stress variation due to exﬁltration is investigated by incorporating the modiﬁed logarithmic bottom boundary layer model of Cheng and Chiew (1998) in the coupled surface-groundwater ﬂow model. The results conﬁrm that as exﬁltration increases, bed shear stress decreases as a result of thickening of the bottom boundary layer.

logarithmic bottom boundary layer model of Cheng and Chiew (1998) in the coupled surface-groundwater flow model. The results confirm that as exfiltration increases, bed shear stress decreases as a result of thickening of the bottom boundary layer.

Introduction
The swash zone is a dynamic region that is alternately covered and uncovered by waves. Coastal processes in the swash zone are particularly difficult to measure due to its frequently turbulent and energetic nature, and not all swash zone processes are fully understood yet. Beach porosity and permeability allow fluid in the coastal aquifer to percolate into or out of the swash zone, exchanging mass and momentum, thereby influencing the swash zone hydrodynamics and morphodynamics and further complicating the dynamics.
As waves run up and down on the beach, the groundwater flow is affected by the wave motion because of the fluctuating water surface (potential head). During the uprush phase, infiltration is the dominant phenomenon (Masselink and Li, 2001;McCall et al., 2014) which has the effect of reducing the momentum of the wave. Greater infiltration rates occur when the water table is lower and the sediment is coarser (Bakhtyar et al., 2011). During the backwash phase, seepage out of the beach is usually observed near the toe of the beach (Li and Barry, 2000;Masselink and Turner, 2012). Seepage from the bed causes modification of the bed shear stress due to its effect on the evolution of the boundary layer. It also influences the effective weight of a sediment particle and hence sediment mobility (Turner and Masselink, 1998;Baldock and Nielsen, 2009;Briganti et al., 2016).
Theoretical analysis and experimental results indicate that exfiltration would cause bed shear stress to decrease due to thickening of the boundary layer, and that infiltration would cause bed shear stress to increase due to thinning of the boundary layer (Cheng and Chiew, 1998;Chen and Chiew, 2004;Corvaro et al., 2014;Miozzi et al., 2015). In addition to in/exfiltration, small-scale flow structures (such as particle inflow/outflow at the seabed) have also been found to be linked to boundary layer dynamics (Corvaro et al., 2014;Miozzi et al., 2015). Cheng and Chiew (1998) and Chen and Chiew (2004) derived modified logarithmic laws for the boundary layer under exfiltration/infiltration, which were validated against experimental results. The numerical study by Pintado-Patiño et al. (2015), which solves the Volume-Averaged Reynolds-Averaged Navier-Stokes equations, further suggests that the effect of seepage on the bottom boundary layer could be so influential that it predominates over bottom boundary layer development, hence the modified logarithmic law developed by Cheng and Chiew (1998) and Chen and Chiew (2004) may not be applicable.
The fluid percolating into or out of the beach can be significant when the swash zone forms part of a larger coastal system such as a barrier system (see Fig. 1). Barrier systems are natural coastal geomorphological features which usually contain porous media connecting the sea to a lagoon (Carter et al., 1989;Carter and Orford, 1993). In a coastal barrier system, different water levels can be present in the beach side and back-barrier lagoon. The water level changes in the lagoon can induce groundwater dynamics near the beachface, which can subsequently affect seepage flow into (or out of) the beach. This exchange of water can subsequently influence hydrodynamical and morphodynamical behaviour at the beach. Therefore, it is important to gain insight into groundwater flow dynamics in a coastal barrier system. Furthermore, such a coastal environment provides a convenient opportunity for the swash zone groundwater flow to be measured in a laboratory setting so that numerical models can be validated.
The prototype-scale BARDEX II (BARrier Dynamics EXperiment II) was conducted in the Delta Flume in the Netherlands in 2012 to understand swash flow, sediment transport and groundwater flow dynamics in a coastal barrier system (Masselink et al., 2016). The data-set acquired in this experiment allows us to test a numerical model in this setting.
Numerical models for both infiltration and exfiltration processes appli-  Figure 1: Schematic diagram of a coastal barrier system, where q is seepage through the beach surface, q x is horizontal pore velocity, q z is vertical pore velocity, u is horizontal depth-averaged water velocity, and n is the normal direction to the boundary surface.
cable in various conditions of beach porosity and permeability have been developed. McCall et al. (2012) and McCall et al. (2014) presented a quasi-3D (depth-averaged) groundwater flow model named XBeach-G for gravel beaches, in which a modified Darcy law was used to describe non-laminar flow due to the high permeability of gravel beaches. Although the model was depth-averaged, the vertical profiles of groundwater head and velocity were estimated following a quasi-3D modelling approach. Additionally, van Gent (1994) developed a coupled surface-subsurface flow model for highly permeable breakwater structures using Forchheimer law. Although the work is mostly related to porous flow in coastal structures, the model could also be applied to gravel barrier beaches. Another coupled surface-subsurface flow model is that by Steenhauer et al. (2012), which modelled air pressure as well as high seepage rates within gravel beaches, using Forchheimer law. In their work, the infiltration/exfiltration at the beach is idealised to be fluid flow driven by a pressure gradient in piston-like movement. For very fine sand (62.5 µm ≤ d 50 ≤ 125 µm, where d 50 is median grain size), Nielsen (1990) and Kang and Nielsen (1997) used the one dimensional Boussinesq equation, the main assumption of which is that the groundwater flow in a shallow aquifer is horizontal, thus neglecting any vertical flow (using the Dupuit-Forchheimer assumption). For beaches that comprise very fine sandy material, which are relatively impermeable, the Boussinesq equation could provide a sufficient description of the flow. However, in porous media of fine to coarse sand (125 µm ≤ d 50 ≤ 1 mm), where flows are laminar and vertical flows are non-negligible, the two dimensional (x − z) Laplace's equation, which is derived based on Darcy law, has been shown to provide reasonable description of the flow (Li et al., 1996(Li et al., , 1997Li and Barry, 2000).
In the present work, we aim to develop a surface-groundwater flow model for a sandy barrier beach. Hence we limit our model to Darcian flow, and solve Laplace's equation to simulate groundwater flow following Li and Barry (2000). One main difference of the present work from Li and Barry (2000) is that a comprehensive bottom boundary layer model (Cheng and Chiew, 1998), which considers the effects of exfiltration, is incorporated in the coupled surface-groundwater flow model, while the Chèzy law was employed in Li and Barry (2000). We use Non-linear Shallow Water Equations (NSWEs) to describe the surface flow, and solve them using a TVD-MacCormack (TVD-MCC) scheme (Briganti et al., 2012). The prototype-scale BARDEX II experimental measurements of surface-groundwater flow are utilised for the model validation. The groundwater flow dynamics and surface-groundwater interactions are examined in detail. We further examine the effects of exfiltration on bottom boundary layer evolution and also on bed shear stress using log-law model in Cheng and Chiew (1998), which is applicable for exfiltration only. As the first step of model development, we assume the beachface is non-erodible to simplify the scenario.
The sections of the paper are organized as follows: the numerical model development is presented in § 2. Then the prototype-scale experiment of BARDEX II as well as the corresponding numerical problem set up are presented in § 3. The model verification and validation are presented in § 4. The effects of exfiltration on the bed shear stress during a single swash event are investigated in § 5. Finally, conclusions are drawn in § 6.

Numerical model
2.1. Groundwater flow model 2.1.1. Governing equations Using the Navier-Stokes equations for mass and momentum conservation, along with Darcy's law for two dimensional groundwater flow, the commonly used Laplace's equation for porous medium flow is obtained (Bear, 2013): where ∇ 2 = ∂ 2 ∂x 2 + ∂ 2 ∂z 2 is the Laplacian operator and φ(x, z, t) is the potential head (m) at horizontal cross-shore position x (m) and vertical position z (m) at time t (s).

Boundary conditions for the groundwater model
For a coastal barrier system, there are generally four boundaries enclosing the groundwater flow domain: AB, BC, CD and DA as shown in Fig. 1.
In the present work, the exit point of the phreatic surface at the sea-side is assumed to be identical to the swash tip following Li and Barry (2000). Therefore, on the left side of the domain (AB), where z b (x) is bed level (m), and h(x, t) is water depth (m). With this approach, transient disconnections between the swash tip and phreatic surface, i.e. seepage faces, will not be captured by the model. However, a similar feature can be captured. The phreatic surface is represented by BC, and the position is described by z = η(x, t). On this surface, pressure p = 0 is assumed, and therefore: The governing equation of the moving phreatic surface is (Liggett and Liu, 1983, see also Appendix B) where p b (dimensionless) is the bed porosity, and k (ms −1 ) is Darcy hydraulic conductivity. We assume that k and p b are both constants. On the face CD, φ = h lag (t).
To test the model performance against theoretical and experimental studies (Kazemzadeh-Parsi and Daneshmand, 2012;Masselink et al., 2016), we assume the barrier beach is founded on an impermeable surface, i.e.: where ∂φ ∂n is the normal derivative of the potential head.

Seepage
With the boundary conditions in § 2.1.2, φ and ∂φ ∂n at the boundaries are either known or computed from the Boundary Integral Equation Method (BIEM) (Liggett and Liu, 1983). This allows the seepage exchange at the beach surface AB and lagoon side CD to be calculated using Darcy's law: where q (ms −1 ) is seepage into (or out of) the boundary (q > 0 indicates exfiltration into the sea or into the lagoon and q < 0 indicates infiltration into the beach). In the case of the lagoon-side boundary, seepage exchange at the permeable CD boundary results in increase or decrease of h lag in time.
After obtaining the boundary solutions, the interior φ, ∂φ ∂x and ∂φ ∂z are calculated using the BIEM scheme. These potential head derivatives allow the pore velocities within the barrier along horizontal and vertical directions (i.e. q x and q z ) to be computed:

Surface flow model
The NSWEs are modified after Wurjanto and Kobayashi (1993) and Clarke et al. (2004) to account for seepage q which results in exchange of mass and momentum between surface and groundwater flow: where τ b is bed shear stress (Nm −2 ), ρ w is water density (kgm −3 ), g is gravitational acceleration (ms −2 ). It is assumed that the surface water density and the groundwater density are the same. Seepage q is calculated from Eq. (7) and the coefficient E in Eq. (10) equals to 1 if q < 0 and 0 if q > 0 because exfiltration is assumed not to contribute to the momentum flux of swash flow (Li et al., 2002).

Bottom shear stress
The shear stress on the beach surface in Eq. (10) is calculated using the momentum integral method in which the horizontal velocity profile inside the BBL is modelled (Clarke et al., 2004;Briganti et al., 2011). We follow Cheng and Chiew (1998) and Clarke et al. (2004) and modify the log-law to include exfiltration. The horizontal velocity within the boundary layer U (x, z , t) on a permeable rough turbulent surface suggested by Cheng and Chiew (1998) is: where z is the vertical distance from the bed surface, u * = |τ b |/ρ w is the friction velocity (ms −1 ), κ is von Karman's constant (we take κ = 0.4) and z 0 is the vertical distance (m) from the bed at which U is assumed to be zero. When there is exfiltration, z 0 is not a constant value, and according to Cheng and Chiew (1998), z 0 is defined as: where K n is bed roughness (m). The variable B in Eq. (12) is a function of q: When there is no exfiltration, q = 0 and B = 8.5. When there is exfiltration, q > 0, and B decreases. Therefore, z 0 increases, which is consistent with the experimental results that the bed is lifted under exfiltration (Cheng and Chiew, 1998;Krogstad and Kourakine, 2000). Eq. (11) is substituted into the BBL momentum equation (Fredsøe and Deigaard, 1993): where U 0 is free stream velocity (derived in Appendix A). Following the approach of Clarke et al. (2004), an ordinary differential equation in Z is formulated: where Z is a function of z 0 and the local boundary layer thickness δ: Then the value of Z is used to determine the friction velocity u * : Finally, the bed shear stress τ b is calculated:

Numerical scheme
The TVD-MCC solver for the NSWEs of Briganti et al. (2012) is used here for the flow above the beach. MCC scheme is a finite difference method where piecewise constant water depth and bed levels are stored at cell centres, and the TVD-function operates as a typical finite volume solver with its nodes coinciding with the cells' centres in the finite difference framework.
The groundwater flow (i.e. Eq. (1)) is solved by the BIEM scheme (Liggett and Liu, 1983). In the BIEM, either φ or ∂φ/∂n needs to be known at each boundary so that the numerical scheme can form a balanced system of equations with equal number of unknown terms (either φ or ∂φ/∂n) and equations so that the unknown term can be computed. At the phreatic surface, both φ and ∂φ/∂n are unknown at a given time interval. Therefore, an additional equation is developed to describe the phreatic surface (Eq. (4)) so that we can have equal number of unknowns and equations. To do this, Eq. (4) is descritized using finite difference and solved along side other equations formed for the other 3 boundaries. For more details regarding the BIEM procedure, the reader is referred to Liggett and Liu (1983). Given by Eq. (3), once φ n+1 is computed, the phreatic surface elevation η n+1 is obtained.
In the BBL sub-model, Z is firstly obtained by solving Eq. (16), which requires U 0 , B and q to be known at t j and t j−1 . U 0 , B and q at t j are from the surface flow model and groundwater flow model. The temporal derivatives of U 0 , B and q are approximated as: dB dt dq dt An adaptive explicit fourth order Runge-Kutta scheme (Press et al., 1989) is used to solve the ODE (15) to get Z at jth time step. As Z is known at each time step, it is used to calculate δ using Eq. (16), u * by Eq. (17) and τ b by Eq. (18). The surface flow model is weakly coupled with the groundwater flow model in the sense that the BIEM is first solved using hydrodynamic data from the previous time step to provide ex/infiltration q, and then hydrodynamic data are updated using the newly calculated q and the hydrodynamic data from the previous time step. However, the NSWEs and Laplace's equation are solved at every time step. We take a constant spatial interval ∆x for TVD-MCC cells, and the BIEM discretization is such that the element size is ∆l cos α = m∆x, where tan α = dz b /dx and m is an integer multiplier of ∆x. Linear interpolation is used to obtain seepage values at every cell for the sea-side boundary (BC in Fig. 1).
The flow chart in Fig. 2 illustrates how these two individual schemes are coupled. Note that in the paper, we refer to water depth in the ith cell at the jth time step as h(i∆x, j∆t) ≡ h j i , with similar expressions for other dependent variables.

Description of the laboratory setup
The BARDEX II experiment was designed to investigate the response of a sandy beach/ barrier system to energetic waves and varying water levels.

BIEM scheme
Runge−Kutta scheme (4th order) The experiment was carried out over 3 months from May to July 2012 at the Delta Flume in the Netherlands (Masselink et al., 2016). The experiment was a follow up to the BARDEX experiment on a gravel barrier beach (Williams et al., 2012). One of the main objectives of this experiment was to investigate the effects of groundwater flow on the swash zone hydrodynamics and morphodynamics in a coastal barrier system. locimeters, 8 pressure transducers (PTs) on the beachface, 15 PTs buried in wells, 4 pairs of PTs in piezometers within the barrier and further 2 PTs in the lagoon region of the coastal barrier system. A complete description of the experimental setup, test conditions and instrumentation was provided by Masselink et al. (2016), here we focus on the instruments used in the current work; hence, only the velocimeters and PTs located in x ≥ 67.5 m are shown in Fig. 3. Irregular random waves generated using a JONSWAP spectrum were produced from a wave maker located at x = 0 m. The significant wave height of the irregular waves was H s = 0.8 m and peak period was T p = 8 s.
Test series A of the BARDEX II experiments is selected for the present work because it was designed to investigate groundwater flow dynamics with different water levels in the lagoon. The focus of this work is on wave run A2 (lagoon level < sea level) and run A4 (lagoon level > sea level) as presented in Table 1 and Fig. 4.

Initial and boundary conditions
The initial sea level and lagoon level for A2 and A4 experiments are illustrated in Fig. 4. The initial free surface elevation at the sea side is set to be 3 m and motionless (u = 0 ms −1 ) before the driving signal is applied from the offshore boundary, located at x = 67.5 m. The initial phreatic surface elevation is unknown, hence it is assumed to be a straight line which connects the sea level and the lagoon level as shown in Fig. 4.
For the runs without waves, the water depth and velocity are kept unchanged at the offshore boundary x = 67.5 m. For the runs with waves, the water depths above the bed were obtained from the calibrated pressure transducer data. Standard calculation methods described in Tucker and Pitt (2001) were used to correct the attenuation of pressure variations with depth and obtain the driving water depth signal at the offshore boundary. The velocity driving signal was obtained by averaging horizontal velocity data from velocimeters elevated at z = 1.83 m, 2.03 m and 2.23 m at the offshore boundary.

Parameter settings
The The beach has porosity, roughness and permeability equivalent to medium sand of d 50 = 0.43 mm. The range of porosity suggested by Masselink et al. (2016) is between 0.37 ≤ p b ≤ 0.42. In the present work, we take p b = 0.4. Following the work of van Rooijen et al. (2012), the bed roughness constant is assumed to be K n = 2.5d 50 = 1.075 × 10 −3 m. The hydraulic conductivity k used is inferred from 20 constant flux tests carried out by Turner et al. (2016) which yielded an average value of 8 × 10 −4 ms −1 . In one wave series of BARDEX II, which lasted nearly 50 mins, the maximum deposition and erosion were in the order of 10 cm. Therefore, the maximum bed change under 300 s considered in the present work is likely to be in the order of a centimeter. Hence, it is reasonable to assume that the bed is non-erodible.

Model verification and validation
4.1. Verification of BIEM scheme: groundwater flow through rectangular dam In this section, the well-known phreatic surface flow case through a rectangular dam (Lacy and Prevost, 1987;Oden and Kikuchi, 1980) is carried out to test the performance of the groundwater flow model with the BIEM scheme: see Fig. 5. Initially, the water levels are fixed on both sides of the dam; the water level on the left side is 10 m whereas the water level on the right is 2 m.
The initial phreatic surface is assumed to be a straight-line that connects the higher water level to the lower water level as shown with a dotted blue line in Fig. 5. An element size of ∆l = 0.1 m was used for the phreatic surface. Both left and right boundaries are Dirichlet boundaries, where φ is known at each node (nodes on the left side φ L = 10 m and the right side φ R = 2 m), and the bottom boundary is of Neumann type, where the flow through the boundary is known to be zero (as ∂φ/∂n = 0). The model is run for 50 s and the phreatic surface is allowed to settle under the difference in potential head between the left and right sides of the dam. The phreatic surface converges to a convex shape as shown in

Validation against the BARDEX II experiments
The test series A2 and A4 of BARDEX II are simulated by the model, and the numerical results are compared against the experimental results.

Groundwater flow without waves
The water table in the model is allowed to settle due to the difference in potential head between the sea side and the lagoon side, without the action of surface water waves, for a period of 300 s. The phreatic surface in the numerical model is allowed to converge for 300 s and comparisons are made against the experimental results in the BARDEX II in Fig. 6 for runs A2 and A4.
The root-mean square deviation between the present model results and the experiment results is computed as: whereη i is the converged/averaged phreatic surface height in the ith node in the model, η i is the corresponding averaged phreatic surface height measured in the experiment and N nodes is the total number of phreatic surface nodes. The comparisons in Fig. 6 show excellent agreement, with computed to be 0.0369 m for case A4 and 0.0324 m for case A2. Contour plots are produced from the converged potential head values and orthogonal lines have been drawn perpendicular to the contours to show the pore water flow in Fig. 6.
As there were no direct measurements of pore velocities below the beach, comparisons are made against inferred time-averaged pore velocities derived from hydraulic head measurements obtained by 4 pairs of PTs located between 83 m < x < 91 m and 0.5 m < z < 2.5 m above the flume floor . The pore velocities are shown at a horizontal spacing of 1 m and a vertical spacing of 0.25 m in Fig. 7, which shows groundwater flow moves from high to low potential head locations. The comparison in Fig. 7 shows good correspondence between numerical and experimental results at x = 89, 90 and 91 m for both A2 and A4 runs. At the same time, the correspondence is also reasonably good at locations z < 1.5 m.
However, discrepancies are observed at x = 86, 87 and 88 m near the beach face where the magnitudes of the velocity vectors are over-predicted in both A2 and A4 runs. These discrepancies suggest that the exfiltration is over-predicted in A2 while infiltration is over-predicted in A4. Near the beach surface, less than 100% saturation (trapped air pockets) could be present and hence this could lead to complex/variable permeability k. In the present simulations, the numerical model assumes a homogeneous bed with a constant permeability which does not take into account such effects. 3e-05 m/s 2 .9 9 5 2 . 9 9 2 . 9 8 2 . 9 6 2 .9 4 2. 92 2.9 2.88

Groundwater flow with waves 4.3.1. Surface water flow comparison
Time series of water depth and velocity from the present model are compared against measurements obtained from the pressure transducers and velocimeters respectively, positioned at x = 72.5 m and 77.5 m in Fig. 8. The water depth comparisons in Fig. 8 show very good phase agreement. Although troughs from model results generally agree well, crests are slightly overpredicted against experimental results. This could also be as a result of the PTs underpredicting the crests due to the non-hydrostatic nature of the flow below the crest (Martins et al., 2017).
The equivalent velocities show less satisfactory agreement with velocimeter data (Fig. 9). Some of the velocity discrepancies could be due to waves that are still dispersive in the shallow water regions. The closer agreement at x = 72.5 m than at 77.5 m could be because x = 72.5 m is closer to the offshore boundary. The velocimeter measurements at x = 77.5 m also suffer from spurious oscillations which suggests air bubbles forming during wave collapse in shallower water depth.
Both comparisons of water depth and velocity against experimental results suggest that the numerical model generally overpredicts the global energy of the waves.

Phreatic surface flow and pore water velocities below the beach
When there are waves, the numerical model is run for a period of 300 s with the initial and boundary conditions in § 3.2. The modelled water table elevations are averaged over the time period and compared against averaged experimental results from BARDEX II in Fig. 10.
When the sea level is lower than the lagoon level ( Fig. 10(a)), the water table position close to the sea side is raised relative to its position without waves ( Fig. 6(a)) and the water table appears to be lowered slightly close to the lagoon as shown in Fig. 10(a). The error = 0.0281 m, which shows very good agreement against experimental results.
When the sea level is higher than the lagoon level, the water table close to the sea side is raised by wave run ups to form a hump-like feature due to the time lag between surface flow and groundwater flow ( Fig. 10(b)). The error = 0.086 m, which suggests reasonably good agreement against experimental results. The greater discrepancy in results closer to the lagoon suggests that the model does not predict very well pore water exit through the sand surface on the lagoon side. Time-averaged interior potential head values predicted by bias and greater magnitude compared with the experimental results. This could be related to the varying permeability discussed in § 4.2.1 as well as rapidly reversing vertical flow near the beach surface. Additionally, the differences could also be related to the differences in modelled hydrodynamics from the surface flow model and the measured hydrodynamics. At the same time, in Fig. 11  for average in the simulation and experimental results can result in certain differences in the pore velocity flow observed below the beach. In addition, the pore velocities are directed towards the sea when the lagoon level < sea level (A4) and waves are present. This is because the run up of waves causes a hump-like feature, the level of which is higher than the lagoon level and also the mean level of the sea. Therefore, water flows from the hump-like feature towards both the lagoon and the sea. This also confirms the findings of Turner et al. (2016) that the surface water wave action on the sea side has a greater influence on the direction of the pore water velocities near the beach than the potential head in the lagoon, ultimately determining whether seepage takes place into (or out of) the beach face.

Flow rates through the barrier
The experimental and numerical average horizontal flow rates through the vertical planes close to the sea side and the lagoon side (locations shown in Fig. 4) are compared in Table 2. The cumulative volume of water passing through each panel is collected over the duration of the simulation, and divided by time to obtain average horizontal flow rates. Table 2: Averaged horizontal flow rates through the sea and lagoon cross-sections shown in Fig. 4 (vertical green lines). Positive flow rates are towards the lagoon and negative flow rates are towards the sea. Experimental results from Turner et al. (2016) are included in brackets to allow comparison.
The average flow rates compare well against experimental results presented in Turner et al. (2016) within a factor of 2-3 difference. Such a difference in the field of experimental hydrogeology is reasonable and confirms the general reliability of numerical model for flux estimates .

Single swash event case
In a prototype-scale laboratory experiment, it is difficult to measure accurately water table fluctuations and seepage flow at the beach during a single swash event. One of the reasons for this is that the seepage and changes in water table could be quite small relative to the resolution and precision of the instruments deployed. However, a numerical model allows us to investigate precisely this. So, we now turn to investigate seepage, water table fluctuations and effects of exfiltration on bed shear stress, during a single swash event.

Seepage at the beach and phreatic surface flow
The seepage flow results are shown in Figs. 12 and 13 for a single swash event. When sea level > lagoon level (A4), the phreatic surface rises along with the swash flow during the uprush phase as seen in Fig. 13(a). Infiltration into the beach is seen in the upper swash region during the uprush. During the backwash phase, a feature similar to a seepage face (Steenhauer et al., 2012) is formed in the upper swash region, because the swash flow retreats towards the sea faster than the groundwater. As the phreatic surface close to the beach has greater potential head than the beach surface, seepage flow takes place from higher potential to lower potential resulting in exfiltration out of the beach. When lagoon level > sea level (A2), there is smaller rise in water table during the uprush phase, as seen in Fig. 12. Similar to A4, infiltration is observed in the upper swash zone but at much smaller rates.
We can also see the seepage pattern below a moving bore (discontinuity) for uprush and backwash phases of a single swash event in Figs. 12 and 13. The seepage patterns show that in front of the moving bore exfiltration (∂φ/∂n < 0) takes place, and behind the bore infiltration (∂φ/∂n > 0) occurs because of the discontinuity in water surface across the bore. The bore seepage velocities are observed in the range of 7.7×10 −4 to 1×10 −3 ms −1 . We observe the seepage across the bore front to be much larger than pore velocities observed in Fig. 11 because the bore front potential gradients ∂φ/∂n are by Packwood and Peregrine (1979) and the numerical results of Li and Barry (2000).

Effects of exfiltration on bed shear stress
Only the effects of exfiltration on bed shear stress are investigated, as the boundary layer model from Cheng and Chiew (1998) is used in the present work and it is applicable for exfiltration only.
The BBL sub-model which incorporates exfiltration is first verified by comparing the simulation results for a single swash event in A4 using the present model with q = 0 and those using the model without exfiltration (Briganti et al., 2011) in Appendix C.
The present BBL model results are validated against the Nielsen et al. (2001) developed an empirical trend line from the Conley and Inman (1994) experiments for the relationship between τ b /τ b0 vs. q/u * 0 (the subscript '0' indicates variables at q = 0). The comparison between the present model results and the empirical trend line for exfiltration ventilation parameters (i.e. 0 < q/u * 0 < 0.025) is shown in Fig. 14. The error between the present model results and the trend line is computed to be = 0.0078 indicating excellent agreement. The decrease in bed shear stress due to ventilation of the boundary layer is well modelled by the present BBL model incorporating exfiltration effects. At the same time, while applicability of Nielsen et al. (2001)'s empirical equation is restricted to 0 < q/u * 0 < 0.025, the present model could be used at higher ventilation parameters. Potentially, this could also improve the accuracy of modelling beach morphodynamics.
In A4, x = 72.5 m is selected for further investigation because we observed in the model validation that exfiltration is generally higher in the lower region of the barrier beach. Therefore, variation in bed shear stress due to exfiltration could also be greater there.
As a bore arrives at x = 72.5 m, the change in ∂φ/∂n across the bore causes alternate exfiltration and infiltration regions in front of and behind the bore respectively (see § 5.1). As the present BBL sub-model incorporates only exfiltration, the seepage transferred from the groundwater model to the BBL is limited to exfiltration only (q > 0) as shown in Fig. 15. Arrival of the bore causes a rapid increase of bed shear stress as the water velocity increases rapidly at that location. At the same time, a rapid increase of exfiltration occurs, but soon after the bore front passes x = 72.5 m shoreward, infiltration takes place behind the bore front as observed in Fig. 13. Infiltration is not transferred to the BBL sub model from the groundwater model as negative q is set to 0. So very rapid increase of dq/dt follows dq/dt → 0 during the period of infiltration. In addition, dB/dt also behaves similarly with B = 8.5 when q = 0. The remaining non-zero terms in Eq. (15) govern τ b and δ for the remainder of the decay observed in Fig. 15 between 16.75 s ≤ t ≤ 18.75 s. Fig. 15 shows that when exfiltration increases, the bed shear stress reduces as a result of thickening of the bottom boundary layer δ. The bed shear stress τ b changes from negative to positive during flow reversal (see Fig. C.17 in Appendix C), because of the change of velocity direction.
Additionally, the BBL model is also compared against Chèzy model, i.e. τ b = −0.5ρ w C d u|u|, using constant friction coefficient C d = 0.01 which is consistent with direct measurements of Barnes et al. (2009) and inferred lab measurements at prototype-scale on rough slopes in Briganti et al. (2011). Both models show that the bed shear stress is maximum at the point of bore arrival, however, the BBL bed shear stress was reduced due to the ventilation of the boundary layer. During the backwash, the Chèzy bed shear stress is greater than BBL shear stress as the Chèzy model assumes a constant friction coefficient C d through the entire duration of the swash event which is greater than the friction coefficient computed from the BBL model.

Conclusions
In the present work, a numerical model that can simulate simultaneously both surface and groundwater flow has been developed and it is validated against results obtained from the prototype-scale BARDEX II laboratory experiment.
Based on the simulation results, the following conclusions can be drawn: • The phreatic surface comparisons against experimental data showed excellent agreement for both A2 and A4 cases without waves with rootmean square deviation 0.03 m ≤ ≤ 0.04 m. The model results with waves showed reasonable agreement against experimental results with 0.03 m ≤ ≤ 0.09 m. The increase in discrepancy when waves are considered arise as a result of differences in free surface depth on the beach, which provide potential head input to the left side boundary of the groundwater model. When the sea level > lagoon level with waves, a hump-like feature forms near the dry region of the beach, under which a flow divide forms. Divergence of pore water from the flow divide towards the sea side and the lagoon side is observed which is consistent with experimental findings in Turner et al. (2016).
• The surface water depth and velocity are well described by the surface water flow model; however, the water depth and velocity results associated with wave crests are not fully satisfactory, but does provide reasonably accurate input to compute potential head within the barrier.
• In the experiments without waves, pore water velocities are always directed from higher potential head to lower potential. This trend reverses when there are waves for the case of sea > lagoon, suggesting that the significant driving action of pore velocities near the beach is the action of surface water waves instead of the difference of sea level and lagoon level.
• During the uprush phase of a single swash event, infiltration is seen to be taking place near the upper swash region. Furthermore, alternate regions of infiltration and exfiltration occur behind and in front of propagating bore respectively, as a result of varying normal potential gradient across the bore front.
• The increase of exfiltration rate leads to the increase in the thickness of the bottom boundary layer, which subsequently results in smaller velocity gradients leading to decrease in bed shear stress. The comparison of the present BBL model shear stress against the Nielsen et al. (2001) empirical trend line showed excellent agreement with an error of = 0.0078.
In the present work, a comprehensive numerical model was presented which can show local swash in/exfiltration reversal from the global effect of the lagoon-sea level changes. The present numerical work was carried out assuming a fixed beach, mainly to allow us to focus on physical processes connected with seepage. Therefore, the seepage effects on the effective weight of sediment particle and critical shield number are avoided as the sediment particle remains fixed on the beach. Future work will include analysis of the effect of seepage on sediment mobility and evolution of the beach morphodynamics. Furthermore, effects of infiltration on the bottom boundary layer also need to be investigated. This will require an alternative bottom boundary layer model as suggested in Chen and Chiew (2004).

Acknowledgement
This work was carried out as part of Eranda Perera's PhD study in the Faculty of Science and Engineering at the University of Nottingham Ningbo China, with funding provided by the University of Nottingham Ningbo China and UK. The authors EP and FZ also would like to acknowledge the financial support from Natural Science Foundation of China (Grant No. 51509135) and Ningbo Science and Technology Bureau, China (Grant No. 2014C50011). The authors also wish to thank Dr. G. Rau for assistance with the interpretation of flow velocity data obtained during the BARDEX II wave flume experiment.
which can be integrated to form: Re-arranging and simplifying the above equation:

Appendix B. Derivation of governing equation for the phreatic surface
Phreatic surface refers to the water surface found within the ground where the gauge pressure is assumed to be zero and a capillary fringe is usually located above it. Here we neglect the presence of moisture above or outside the free surface (i.e. capillary fringe), and assume that the free surface is an abrupt interface between air and water in the void space.
For saturated groundwater flow within a barrier, the phreatic surface elevation can be represented by z = η(x, t). As the pressure at all points on the phreatic surface is taken as p = 0, from φ = z + p/(ρ w g): φ = η(x, t) on z = η. (B.1) The unsteady (moving) phreatic surface is a fluid surface which always contains the same fluid particles (Bear, 2013). Therefore, the rate of change of the position of the phreatic surface must be equal to the vertical velocity on the surface: where q/p b is the velocity of the fluid particles on the phreatic surface and p b is the void space available in unit area of porous medium which allows water to percolate. The present BBL model is compared with numerical results from the BBL solver in Briganti et al. (2011) which does not include seepage effects. The verification is performed by running the present model by setting permeability k = 0 ms −1 (which results in q = 0 ms −1 ). When q = 0 ms −1 , dB dt = dq dt = 0 in Eq. (15), which reduces to: Eq. (C.1) is the same ordinary differential equation solved in Briganti et al. (2011). The difference between the models lies in the definition of z 0 which is assumed to be z 0 = K n /30 instead of Eq.