Hydrodynamic Modeling of a Reef-Fringed Pocket Beach Using a Phase-Resolved Non-Hydrostatic Model

The non-hydrostatic wave-flow model SWASH was used to investigate the hydrodynamic processes at a reef fringed pocket beach in southwestern Australia (Gnarabup Beach). Gnarabup Beach is a ~1.5 km long beach with highly variable bathymetry that is bounded by rocky headlands. The site is also exposed to large waves from the Southern Ocean. The model performance was evaluated using observations collected during a field program measuring waves, currents and water levels between June and July 2017. Modeled sea-swell wave heights (periods 5–25 s), infragravity wave heights (periods 25–600 s), and wave-induced setup exhibited moderate to good agreement with the observations throughout the model domain. The mean currents, which were highly-spatially variable across the study site, were less accurately predicted at most sites. Model agreement with the observations tended to be the worst in the areas with the most uncertain bathymetry (i.e., areas where high resolution survey data was not available). The nearshore sea-swell wave heights, infragravity wave heights and setup were strongly modulated by the offshore waves. The headlands and offshore reefs also had a strong impact on the hydrodynamics within the lagoon (bordered by the reefs) by dissipating much of the offshore sea-swell wave energy and modifying the pattern of the nearshore flows (magnitude and direction). Wave breaking on the reef platforms drove strong onshore directed mean currents over the reefs, resulting in off-shore flow through channels between the reefs and headlands where water exchanges from the lagoon to ocean. Our results demonstrate that the SWASH model is able to produce realistic predictions of the hydrodynamic processes within bathymetrically-complex nearshore systems.


Introduction
Reefs (including coral and rocky) are common features in nearshore environments worldwide. The physical mechanism underlying the hydrodynamic process on reef beaches can be similar to sandy beaches, however the presence of the reefs introduces key differences. For example, infragravity (IG) waves (with periods of 25-600 s) have been shown to often account for a considerable portion of the total wave energy, e.g., [1,2] and thus are key to developing an understanding of the hydrodynamic processes. Reefs are often also detached from the shoreline and/or are discontinuous in the alongshore thus introducing two-dimensionality in the waves and the resulting wave-driven mean currents. This two-dimensionality is often accompanied by large bottom roughness and steep slopes formed by the reef morphology. While The model is then used to understand the wave and circulation dynamics along the beach during different wave and water level conditions.

Study Site and Field Observations
The study focused on Gnarabup Beach in the southwest of Western Australia (Figure 1) [29]. The beach is a 1.5 km long pocket beach that is fronted by two limestone reefs located~600 m offshore. The reefs become partially exposed at low water levels (see Figure 1a), forming a semi-protected lagoon (hereinafter referred to as 'lagoon' for simplicity). The area behind the reefs has depths ranging from 2 to 5 m with the seabed consisting of a mix of medium to coarse carbonate sands, seagrass, and scattered limestone rocks. With exposure to the Southern Ocean, the offshore wave conditions consist of mostly swell, with offshore significant wave heights averaging~2-3 m but can reach~8 m during storm conditions [29]. The tides at the site are micro-tidal with average and maximum tidal range during the study period of 0.36 m and~1.25 m, respectively.
A 4-week field study was conducted in June-July 2017 to quantify the transformation of waves from offshore to the beach and to record the circulation dynamics within the lagoon as reported in [29]. During the study, an array of 22 instruments were deployed offshore and inside the lagoon (Figure 1a). Two Nortek AWACs were placed outside the lagoon at~19 and~21 m depth directly offshore of each reef. Inside the reef system, 8 current meters (4 Nortek Aquadopp Profilers, 2 Nortek Vectors, and 2 Teledyne RDI Acoustic Doppler Current Profilers (ADCP) and 12 RBRsolo pressure sensors were situated in cross shore and alongshore arrays to record the waves, water levels and currents ( Figure 1a, Table 1).
Directional wave spectra from the offshore AWAC sensors (S0 and C0) were computed with the Maximum Likelihood Method using the Nortek Storm software. Inside the lagoon, wave spectra (converting from the pressure spectra using linear wave theory) were calculated at each sensor location from hourly records of pressure collected at 1 or 2 Hz (3600 or 7200 samples, see Table 1) using Welch's averaged periodogram method with a Hamming window of 50% overlap to reduce spectral leakage. The wave spectra were integrated within the sea-swell band with periods of 5-25 s (0.2-0.04 Hz) and infragravity (IG) with periods 25-600 s (0.0017-0.04 Hz) to compute significant wave heights within each band. The observed wave setup for the lagoon sensors was estimated from depth difference between the sensors located at the reef and offshore [29][30][31], where η is the wave setup at the lagoon sites; h and h 0 are the mean depth obtained from the recorded pressure time series at the lagoon site and at the offshore AWAC-C0 respectively, with the overbars indicating hourly averaging; and ∆h is the difference in elevation between the offshore sensor and the lagoon sensors. In this study, we neglected setup calculations at sites C1, C4 and S7, due to evidence of some drift in the setup estimates, possibly caused by either sensor drift or the instruments settling during the deployment. The velocity profiles recorded by the acoustic current profilers were depth-averaged and then time-averaged (hourly). Also, we assumed the velocity information from the Nortek Vectors, which provided point velocity measurements just above the sensor located~1 m above the bottom, to be comparable to depth averaged currents due to the relatively shallow depths (2.5-3 m, see Table 1). In this study, we focused on swell events (T p >~10 s, see Figure 2a     Cyclic boundary conditions were applied along the lateral boundaries and were placed ~700 m away from the headlands that enclosed the pocket beach (see Figure 1c). With this approach, the model domain was continuous in the alongshore direction. To ensure that the bathymetry is repetitive without alongshore discontinuities, the domain was extended with a 300 m gradual transition between the northern and southern lateral boundary. The final model domain was 1.8 km by 3.0 km, resulting in spatial grid of 600 × 762 points with minimum and maximum bottom

Numerical Model Description
SWASH numerically solves the Reynolds-averaged Navier-Stokes (RANS) equations for an incompressible fluid using second order finite difference scheme see [18], and references therein for details regarding the numerical implementation. For a three-dimensional flow field in the x (cross-shore), y (alongshore), z (vertical) directions, the governing equations of the model are: t is time, ζ is surface elevation measured from still water level d, h = ζ + d is the total water depth, u, v and w are the respective flow velocities in x-, yand zdirection (with uppercase indicating depth-averaged velocities), g is the acceleration due to gravity, q is the non-hydrostatic pressure, and τ represents the turbulent stress terms. At the bottom boundary, a quadratic friction law was included to model the stress induced by the seabed on the flow, in which c f is a dimensionless friction coefficient. For the turbulent stresses, separate eddy viscosities are used for horizontal and vertical mixing see [32] for details. The horizontal eddy viscosities in this work were computed from a Smagorisnky type formulation (with default parameters), and a constant vertical viscosity was used to enhance vertical mixing (1 × 10 −5 m 2 /s).

Model Set-Up
SWASH was used to simulate the hydrodynamic conditions at Gnarabup Beach that occurred between 30 June and 3 July 2017 during 17 distinct hourly periods that were representative of the range of conditions. This approach was used to balance capturing a wide range of conditions with the highly-resolved (temporally and spatially) simulations over the large domain; continuous simulations over the entire period were not computationally feasible even using the high performance computing resources available. The numerical domain for Gnarabup Beach was constructed using a uniform rectangular computational grid with 3 m resolution in the cross-shore and 4 m in the alongshore, which was rotated 11.7 • clockwise to align with the dominant direction of the coastline. This grid size was determined through a sensitivity analysis to grid resolution described in Appendix A. Bathymetry for the area was obtained through the Western Australia Department of Transport, with a majority of the area covered by a multi-beam survey collected in 2016 and some areas only covered by single beam survey from 1995 (refer to the spatial coverage of the available bathymetry in Figure 1b). Beach topography data were obtained from a backpack-based differential Global Positioning System (GPS) survey collected in February 2017 (see Figure 1b). While there was generally good bathymetry coverage over much of the study area, bathymetry data was not available for the very shallow areas over the reefs and headlands.
The depths in these areas were estimated at 0 m (relative to the Australian Height Datum, AHD, approximately mean sea level) for the areas that appeared exposed on aerial imagery collected during low water conditions. Adjacent areas were approximated with depths of −1 or −0.5 m based on the estimated submergence from the aerial imagery. The surveyed bathymetry and topography were linearly interpolated onto a 1 × 1 m grid, which was then used to generate the bathymetry within the model domain ( Figure 1c). The offshore open boundary was located~500 m offshore from the AWAC position, which includes a region with constant water depth of 25 m AHD (to ensure stable wave generation) and a gradual transition to the observed bathymetry.
Cyclic boundary conditions were applied along the lateral boundaries and were placed~700 m away from the headlands that enclosed the pocket beach (see Figure 1c). With this approach, the model domain was continuous in the alongshore direction. To ensure that the bathymetry is repetitive without alongshore discontinuities, the domain was extended with a 300 m gradual transition between the northern and southern lateral boundary. The final model domain was 1.8 km by 3.0 km, resulting in spatial grid of 600 × 762 points with minimum and maximum bottom elevations of −25 m and 6.5 m, respectively.
At the offshore boundary, the 17 hourly scenarios were simulated that are part of a four day period (30 June to 3 July 2017) that include a wide range of offshore wave conditions (H s from~2-5 m) and water levels (Figure 2a-d). In this study we used hourly directional wave spectra estimated from the AWAC at C0, which was imposed uniformly along the offshore open boundary. The wave spectra measured at AWAC C0 were back-refracted and de-shoaled to the location of the wavemaker following the approach of [33]. Based on the directional wave spectra, short-crested sea states were generated in the model using a weakly reflective wavemaker including a second order correction to include bound infragravity waves [25].
The initial time step (∆t) within the computation was set to 0.025 s to satisfy 0.3 > Courant-Friedrichs-Lewy (CFL) < 0.8, with SWASH internally modifying the time step thereafter. Two equidistant vertical sigma (terrain-following) layers were used (i.e., each layer occupied 50% of the local water depth), which is sufficient to capture the dispersive characteristics of the considered sea states. We used the hydrostatic front approximation to capture the initiation of wave breaking [34] given the coarse vertical resolution.
For each of the 17 wave and water level conditions considered (Figure 2a-d), the model was run for 60 min plus an additional 30 min to allow for model spin-up. The spin-up period was defined based on examining the evolution of the integrated kinetic energy (KE), potential energy (PE) and enstrophy (Z, i.e., the integral of vorticity-Ω squared) over the model domain where L y and L x are the cross-shore and long-shore lengths of model domain, respectively ( Figure 3). The model was assumed to be in 'steady state' when these variables fluctuated about a steady mean see [33], for detailed explanation. Output variables consisted of surface water level elevation and depth-averaged velocities (U and V), which were saved at 2 Hz over a 4 × 4 m grid. To minimize the overall computational time (including output generation), variables were only saved for an area of 1 km (cross-shore) × 1.8 km (alongshore) covering the area of interest (Figure 1c). Time series output were also saved at the location of each instrument within the domain (i.e., 22 instrument locations). The SS-IG waves and mean currents from the model were calculated using the same methods applied to the field data in Section 2.1. An example of the model output (instantaneous free surface elevations, SS and IG wave heights, and mean currents) are shown for one of the scenarios in Figure 4, with the estimated setup in Figure 5a.
where Ly and Lx are the cross-shore and long-shore lengths of model domain, respectively ( Figure 3). The model was assumed to be in 'steady state' when these variables fluctuated about a steady mean see [33], for detailed explanation. Output variables consisted of surface water level elevation and depth-averaged velocities (U and V), which were saved at 2 Hz over a 4 × 4 m grid. To minimize the overall computational time (including output generation), variables were only saved for an area of 1 km (cross-shore) × 1.8 km (alongshore) covering the area of interest ( Figure 1c). Time series output were also saved at the location of each instrument within the domain (i.e., 22 instrument locations). The SS-IG waves and mean currents from the model were calculated using the same methods applied to the field data in Section 2.1. An example of the model output (instantaneous free surface elevations, SS and IG wave heights, and mean currents) are shown for one of the scenarios in Figure 4, with the estimated setup in Figure 5a.
Considering the spatial scale of the model, all simulations in this study were run on a supercomputer at the Pawsey Supercomputing Centre (Western Australia). Each one-hour simulation took about 20 h to complete using 48 cores.     Considering the spatial scale of the model, all simulations in this study were run on a supercomputer at the Pawsey Supercomputing Centre (Western Australia). Each one-hour simulation took about 20 h to complete using 48 cores.

Model Performance
Model performance was assessed against the observations at each of the 22 instrument sites using the Willmott Skill (WS), e.g., [35][36][37] which is computed as where X mod is a model prediction, X obs is an observation, and overbars indicate averaging over the 17 scenarios. WS was computed at each site for various bulk wave and flow parameters for all 17 hourly scenarios (N = 17), with the over-bar indicating averaging over the 17 scenarios. Good model performance is considered when WS > 0.6, moderate performance when WS is 0.3 to 0.6, and poor when WS < 0.3. In addition to the WS, we also computed the root mean square error (RMSE) and the mean Bias to assess the model performance,

Model Calibration
The model was calibrated for a single sea-state that represented relatively moderate conditions during the considered study period (H s0 = 3.14 m with a T p0 = 13.4 s at the offshore boundary).
The calibration process focused on assessing the performance to variations in the empirical breaking parameter α [34] and bottom friction coefficient c f based on the comparison of the modeled and observed waves, mean water levels and currents at the 22 instrument sites. Changing the breaking parameter, α, over a range from 0.04 to 3 did not impact the model skill significantly (not shown), so we used the default value of α = 0.6 within SWASH. We next focused on the sensitivity of the model output to the uniform bottom friction coefficient (c f ) to reproduce the hydrodynamic parameters within the model domain. While other drag formulations are available in SWASH (e.g., Manning), we adopted a constant friction coefficient that is commonly applied, including in reef environments, e.g., [27,38]. We varied the c f over a range from 0.01 to 0.08 across the domain in intervals of 0.01 for the first eight simulations and higher intervals 0.02 from 0.08 to 0.16 and 0.05 from 0.16 to 0.3. As we describe in Section 4, we also attempted a spatially-varying c f but did not find improved results for the c f values considered.
A summary of the RMSE and mean Bias between observation and model for SS and IG wave heights, setup η, and currents (magnitude, U and V components) for various bottom friction coefficients is shown in Figure 6, in which the sensors were grouped according to their location within the study site (offshore, channel and lagoon, Figure 1a). For each area, we present the average RMSE and Bias as well as the standard deviations across the sites within each grouping (with the exception of the offshore sites as there were only 2 sensors). At the offshore sites, the error and Bias (blue lines in Figure 6) were small (consistent with these sites being used to derive the boundary forcing). For shallower depths, i.e., in the channels (green lines) and the lagoon (red lines), as c f initially increased the results in general improved for wave heights, setup and currents, reached an optimum level of performance, and then often deteriorated for higher values of c f . Across the range of c f values considered, the overall model errors across all variables (SS and IG wave heights, setup, and currents) were lowest using a c f = 0.05 ( Figure 6). The estimated c f value is within the range of other studies comparing from reef sites typically of the order 0.01-0.1, e.g., [9,[39][40][41]. Thus, the c f value of 0.05 was subsequently applied to the full suite of scenarios.

Model Application
Following selection of an optimal bottom friction coefficient (cf) based on a typical (moderate) wave condition, we evaluated the model for all 17 scenarios over the 4-day period. Below, the

Model Application
Following selection of an optimal bottom friction coefficient (c f ) based on a typical (moderate) wave condition, we evaluated the model for all 17 scenarios over the 4-day period. Below, the predicted SS and IG wave heights, current velocities (U and V vectors, as well as the magnitude) and wave setup are compared against the observations. Similar to the analysis in the previous section, the results were grouped into 3 different zones of offshore, channels and lagoon. We also divided the lagoon into northern, central and southern areas based on the position of each instrument within the lagoon relative to the reefs (see yellow dashed lines in Figure 1a). The northern portion of the lagoon is deeper and relatively exposed to offshore waves (compared to the other parts) with scattered rocks near the shoreline (Figure 4). The central area is located behind two offshore reefs and influenced by strong currents due to wave breaking on the reefs. The southern area is the most protected part of the lagoon and typically has low SS wave energy and weak currents due to protection from the rocky headland.

Wave Height and Setup
Across the 17 scenarios, SWASH predictions showed the typical wave transformation processes within a coastal region. SS waves shoaled and subsequently broke on the reefs (Figure 4a,b). The dissipation of SS waves at the fringing reefs coincided with an increase in the IG wave heights (Figure 4c). Inside the lagoon, SS wave heights were significantly smaller than offshore, and the ratio between IG and SS waves was high (often greater than 0.6), especially near the shoreline, indicating the importance of IG waves in dictating shoreline water levels (not shown). Within the surf-zone (primarily over shallow reef areas), SWASH simulated breaking on the reefs and refraction/diffraction in their lee (Figure 4a). Mean water levels are lower in the offshore region where wave shoaling occurs (i.e., causing set down), and elevated due to wave setup on the reefs and near the shoreline ( Figure 5) coinciding with regions of strong wave breaking.
The performance of SWASH in predicting the observed waves and currents within each area of the study site are summarized in Figures 7 and 8, and Table 2. The model performed well at reproducing the SS wave heights at the offshore site (WS ≥ 0.6, Bias < 5 cm) ( Table 2). In the channels (Figure 7b), we noticed that for most of the 17 scenarios the model over predicted the SS wave heights (e.g., at site C1 see the blue dashed line and green dots in Figure 8a) and S1 (not shown), resulting in moderate skill at Central channel (WS = 0.55) but good skill at Southern channel (WS of 0.64) with a Bias at both locations of~0.5 m. The SS wave heights shoreward of the southern channel (S2 and S3) tended to be better predicted (WS ≥ 0.6, Bias <~5 cm) than those shoreward of the central channel (C2 and C3, with WS~0.3 and Bias > 0.2 m). Within the lagoon, SWASH predicted the SS wave heights in the area behind the southern offshore reefs with good accuracy (C8 and almost all sensors at the south part) with WS fluctuating from 0.64-0.85 and Biases ranging from 0.01-0.25 m. The individual lagoon sites with poor skill (WS below 0.3) mostly occurred at sites behind the northern reefs and the adjacent area (i.e., northern (N3) and central (C5-C6) parts), except at C4 that was located near the northern fringing reef where WS = 0.85. We note that regions with poor SS wave height model performance generally tended to be related to regions where there was the most uncertainty in the bathymetry; for example, in the northern portion of the lagoon the available bathymetry data was limited (see Figure 1b).
Scatter plots between the observed and modeled IG wave heights (Figure 7d-f) showed good agreement between the IG wave heights with WS ranging from 0.62 to 0.9 at all of the sites, with maximum Bias of~0.1 m ( Table 2). The RMSE of the modeled IG wave heights across three areas (northern, central and southern) were relatively small, ranging from 0.03 to 0.12 m for all sensors. At deeper areas, i.e., the AWAC sensors (C0 and S0), the predicted IG waves tended to be less accurate (over-predicted) during larger wave events (Figures 7d,e and 8e), causing the overall WS score at these sensors to be around 0.6-0.7. Inside the lagoon, the IG wave predictions at all sensors showed WS ≥ 0.6 with Biases < 0.05 m. By comparing the IG waves at all sensors, we found that the IG waves showed less spatial variability across the sites (Figures 7f and 8f-h), which indicates the IG waves were less impacted by the variable water depths in the lagoon.
Setup was, in general, well predicted with good WS at most sensors with slight over-prediction by the model compared to the observations, with a mean Bias of 0.05 m across all sites (Table 2). Despite the slight over-prediction of setup at all locations (Figure 7g,h), the model showed fairly consistent Bias and tracked the variations in the observed setup (Figure 8i-k). Further comparison between observation and model for each sensors showed that predicted setup in the northern and central areas (which are deeper and more exposed to offshore waves) was more accurate (WS > 0.8) than that in the more protected southern. This is in contrast with the SS wave heights for which the model performed better (on average) in the southern portion of the lagoon (Table 2).
reproducing the SS wave heights at the offshore site (WS ≥ 0.6, Bias < 5 cm) ( Table 2). In the channels (Figure 7b), we noticed that for most of the 17 scenarios the model over predicted the SS wave heights (e.g., at site C1 see the blue dashed line and green dots in Figure 8a) and S1 (not shown), resulting in moderate skill at Central channel (WS = 0.55) but good skill at Southern channel (WS of 0.64) with a Bias at both locations of ~0.5 m. The SS wave heights shoreward of the southern channel (S2 and S3) tended to be better predicted (WS ≥ 0.6, Bias < ~5 cm) than those shoreward of the central channel (C2 and C3, with WS~0.3 and Bias > 0.2 m). Within the lagoon, SWASH predicted the SS wave heights in the area behind the southern offshore reefs with good accuracy (C8 and almost all sensors at the south part) with WS fluctuating from 0.64-0.85 and Biases ranging from 0.01-0.25 m. The individual lagoon sites with poor skill (WS below 0.3) mostly occurred at sites behind the northern reefs and the adjacent area (i.e., northern (N3) and central (C5-C6) parts), except at C4 that was located near the northern fringing reef where WS = 0.85. We note that regions with poor SS wave height model performance generally tended to be related to regions where there was the most uncertainty in the bathymetry; for example, in the northern portion of the lagoon the available bathymetry data was limited (see Figure  1b).   Figure 1a). Numbers on the subplots represent average values of the respective Willmott Skills, RMSEs and Bias (refer to Table 2). deeper areas, i.e., the AWAC sensors (C0 and S0), the predicted IG waves tended to be less accurate (over-predicted) during larger wave events (Figures 7d,e and 8e), causing the overall WS score at these sensors to be around 0.6-0.7. Inside the lagoon, the IG wave predictions at all sensors showed WS ≥ 0.6 with Biases < 0.05 m. By comparing the IG waves at all sensors, we found that the IG waves showed less spatial variability across the sites (Figures 7f and 8f-h), which indicates the IG waves were less impacted by the variable water depths in the lagoon. Setup was, in general, well predicted with good WS at most sensors with slight over-prediction by the model compared to the observations, with a mean Bias of 0.05 m across all sites (Table 2). Despite the slight over-prediction of setup at all locations (Figure 7g,h), the model showed fairly consistent Bias and tracked the variations in the observed setup (Figure 8i-k). Further comparison between observation and model for each sensors showed that predicted setup in the northern and central areas (which are deeper and more exposed to offshore waves) was more accurate (WS > 0.8) than that in the more protected southern. This is in contrast with the SS wave heights for which the model performed better (on average) in the southern portion of the lagoon ( Table 2).

Depth-Averaged Currents
We compared the modeled mean (hourly) currents from the 17 scenarios with the coinciding observations from the 10 sites that measured velocity located offshore (C0 and S0), in the channels (C1 and S1), as well as central and southern parts of the lagoon (C4, C6, S6, S7, S9, and S10) (Figure 1a). The observed currents across the 17 wave conditions showed a broad variability of flows within the study area, ranging from very weak currents inside the lagoon at the southern corner of the beach (S10) with average of 0.02 m/s (maximum of 0.03 m/s) to strong currents recorded in the channel sites (e.g., with average of 0.74 m/s, maximum of 1.02 m/s, at S1). Evaluation of the mean currents indicated that SWASH tended to under-predict the magnitudes (Figure 9a), except at C6 and S10 where the model over-predicted the current magnitudes. At all sites, we found that the Bias in the current magnitude was less than 0.1 m/s ( Table 2). The model skills for each locations were moderate (0.3-0.6) for predicting currents at most sites, except at site C1 and S1 where there was good skill (WS ≥ 0.6) and one sites with poor skill at S7, as shown in Table 2.
shows that the model generally under-predicted the currents (Figure 10a,c). At the channel sensors, the model captured increasing current strengths for more energetic waves inside the channels (Figure  10a), although current directions were off by about 25° (Figure 10b). The current predictions were found to be in better agreement for larger off-shore wave conditions, with the incident wave directions were relatively constant across all wave conditions (the mean incident wave directions was 251 ± 3°, see Figure 2c).  The model also demonstrated moderate (at C0, S0, S6, and S10) and good skill (C6) in predicting current directions, even though the other sites exhibited poor skill ( Table 2). Note that at some points, the discrepancies between observed and predicted current directions could be very large due to phase wrapping (e.g., >300 • ) resulting in poor skill, even though the comparison (Figure 9a) showed the directions did not differ too much. Hence, for these conditions, we considered the adjustment in direction by 360 • to account for phase wrapping. At the channel sensors (C1 and S1, with poor WS at both sites), the modelled outflow had a more southerly direction compared to the observations. At the offshore sensors the model did not capture the observed northern flows, which is likely explained by larger scale flows that are not captured within the model. Despite these discrepancies, the model captured the overall flow patterns within the lagoon system, with weak flows in the lagoon and strong outflow through the channels (Figure 9a). Figure 10 shows the comparison between the observed and modeled currents at the channels (C1 and S1) and inside the lagoon (C4, C6, S6, S7, S9, S10) for each of the 17 scenarios. This comparison shows that the model generally under-predicted the currents (Figure 10a,c). At the channel sensors, the model captured increasing current strengths for more energetic waves inside the channels (Figure 10a), although current directions were off by about 25 • (Figure 10b). The current predictions were found to be in better agreement for larger off-shore wave conditions, with the incident wave directions were relatively constant across all wave conditions (the mean incident wave directions was 251 ± 3 • , see Figure 2c).

General Circulation Dynamics
Despite the wave heights varying significantly over the study period (from 2.4 m to nearly 5 m), the observed wave directions were relatively consistent (247 to 255°; see Figure 2c) due to the persistent arrival of Southern Ocean swell along this portion of coastline. While the flow magnitudes varied, we observed similar circulation patterns over the 17 scenarios (Figure 9b). From the model results, we observed that wave breaking over the two reefs resulted in large setup on the reefs (up to 0.46 m) which drove onshore flows into the lagoon where setup was lower (Figure 5b). Towards the outer edges of the reefs, the modelled currents diverged directly into the channels resulting in lower mass flux into the shallower portions of the lagoon. Modeled currents along the shoreline (up to ~200 m offshore) were generally weak (averaging <0.06 m/s) in the central and southern portions of the lagoon as a result of the divergence of the flow from the reefs directly into the channels. In the deeper northern portion of the lagoon, the modeled currents were overall stronger along the beach, likely due to a combination of higher wave penetration into the lagoon, the larger gap between the northern reef and northern headland, as well as the increased depths compared to that of the central and the southern area (Figure 9b). Overall, the model revealed that current magnitudes near the shoreline were the largest in the northern portion of the lagoon and weakest in the central portion of the lagoon (Figure 9). The model also indicated some drainage of water from the central portion to the northern portion of the lagoon through the gap between the northern reef and shoreline (Figure 9b), associated with an alongshore pressure gradient (lower setup in northern part lagoon, Figure 5b). In the channels, the flow out of the lagoon was approximately proportional to the offshore wave height indicating the larger onshore mass flux into the lagoon and channels during larger wave conditions ( Figure 11).

General Circulation Dynamics
Despite the wave heights varying significantly over the study period (from 2.4 m to nearly 5 m), the observed wave directions were relatively consistent (247 to 255 • ; see Figure 2c) due to the persistent arrival of Southern Ocean swell along this portion of coastline. While the flow magnitudes varied, we observed similar circulation patterns over the 17 scenarios (Figure 9b). From the model results, we observed that wave breaking over the two reefs resulted in large setup on the reefs (up to 0.46 m) which drove onshore flows into the lagoon where setup was lower (Figure 5b). Towards the outer edges of the reefs, the modelled currents diverged directly into the channels resulting in lower mass flux into the shallower portions of the lagoon. Modeled currents along the shoreline (up to~200 m offshore) were generally weak (averaging <0.06 m/s) in the central and southern portions of the lagoon as a result of the divergence of the flow from the reefs directly into the channels. In the deeper northern portion of the lagoon, the modeled currents were overall stronger along the beach, likely due to a combination of higher wave penetration into the lagoon, the larger gap between the northern reef and northern headland, as well as the increased depths compared to that of the central and the southern area (Figure 9b). Overall, the model revealed that current magnitudes near the shoreline were the largest in the northern portion of the lagoon and weakest in the central portion of the lagoon (Figure 9). The model also indicated some drainage of water from the central portion to the northern portion of the lagoon through the gap between the northern reef and shoreline (Figure 9b), associated with an alongshore pressure gradient (lower setup in northern part lagoon, Figure 5b). In the channels, the flow out of the lagoon was approximately proportional to the offshore wave height indicating the larger onshore mass flux into the lagoon and channels during larger wave conditions (Figure 11).
Despite the reefs not offering complete alongshore coverage (only occupying approximately one-third of the distance between the headlands) they still provide effective protection to the beach areas, especially in the southern and central portions of the lagoon. The model and observations show that the IG wave heights were more uniform across the lagoon (Figure 7f and [29]).
Despite the reefs not offering complete alongshore coverage (only occupying approximately one-third of the distance between the headlands) they still provide effective protection to the beach areas, especially in the southern and central portions of the lagoon. The model and observations show that the IG wave heights were more uniform across the lagoon (Figure 7f and [29]). Figure 11. Modeled mean (hourly) current velocity at South (S1), Central (C1) and North (N1) channels (refer to Figure 1a for the locations) versus observed wave height at the AWAC (C0).

Discussion
In this study, we assessed the ability of SWASH to simulate the hydrodynamics of a complex nearshore site, with rocky reefs fronting a pocket beach bounded by rocky headlands. The model performance varied depending on the hydrodynamic processes that were simulated and spatial location within the study site. SWASH accurately predicted the IG waves (good skill, low RMSE, low Biases) at all the observed locations for 17 (hour) scenarios. It was found to reproduce the SS waves and setup well at majority of the sensors (see Table 2). When predicting the mean wave-driven currents, SWASH exhibited a wider range of skill (poor to good, WS 0.29 to 0.66), with relatively small RMS error (average of 0.09 m/s) and Bias (average of 0.07 m/s).
In our study, we identified some possible factors that potentially influenced the accuracy of the model. The specific locations where the SS waves and setup were most poorly predicted were mainly located on the northern part of the beach (see Table 2). The certainty of the bathymetry in the northern part of the beach is lower due to less survey coverage (see Figure 1b), with the depth in some areas (the shallow reefs) estimated from aerial images. In addition, this region features a number of small patch reefs (finer than grid cells) that are scattered within the domain (mainly at the northern part) could not be fully resolved in the model. These small reefs may cause SS waves to break in different areas than predicted in the model, thus likely contributing to the magnitude and direction of the resulting wave-driven flows, e.g., [42].
Other possible factors that potentially control the accuracy of modeled SS waves and setup is the type of computational grid used in the model. In this study, we used a structured rectilinear grid with a 3 (cross-) × 4 (along shore) m grid resolution. The main drawback of such a structured grid with a relatively coarse grid resolution is the limited flexibility to capture multiscale features of complex bathymetry [43]. As a result, the model likely did not capture all the details of the complex bathymetry of Gnarabup Beach. Furthermore, phase-resolving models like SWASH require fine resolutions to capture higher frequency waves, e.g., [1]. However, the use of finer grid sizes in this study, which greatly increases computational cost, did not significantly improve the results, as shown in Appendix A. This may reflect that to gain the benefits of higher resolution computations requires having sufficiently high-resolution bathymetry at these fine-scales, which was not available in this study. Presently SWASH also only allows uniform wave conditions (based on single directional Figure 11. Modeled mean (hourly) current velocity at South (S1), Central (C1) and North (N1) channels (refer to Figure 1a for the locations) versus observed wave height at the AWAC (C0).

Discussion
In this study, we assessed the ability of SWASH to simulate the hydrodynamics of a complex nearshore site, with rocky reefs fronting a pocket beach bounded by rocky headlands. The model performance varied depending on the hydrodynamic processes that were simulated and spatial location within the study site. SWASH accurately predicted the IG waves (good skill, low RMSE, low Biases) at all the observed locations for 17 (hour) scenarios. It was found to reproduce the SS waves and setup well at majority of the sensors (see Table 2). When predicting the mean wave-driven currents, SWASH exhibited a wider range of skill (poor to good, WS 0.29 to 0.66), with relatively small RMS error (average of 0.09 m/s) and Bias (average of 0.07 m/s).
In our study, we identified some possible factors that potentially influenced the accuracy of the model. The specific locations where the SS waves and setup were most poorly predicted were mainly located on the northern part of the beach (see Table 2). The certainty of the bathymetry in the northern part of the beach is lower due to less survey coverage (see Figure 1b), with the depth in some areas (the shallow reefs) estimated from aerial images. In addition, this region features a number of small patch reefs (finer than grid cells) that are scattered within the domain (mainly at the northern part) could not be fully resolved in the model. These small reefs may cause SS waves to break in different areas than predicted in the model, thus likely contributing to the magnitude and direction of the resulting wave-driven flows, e.g., [42].
Other possible factors that potentially control the accuracy of modeled SS waves and setup is the type of computational grid used in the model. In this study, we used a structured rectilinear grid with a 3 (cross-) × 4 (along shore) m grid resolution. The main drawback of such a structured grid with a relatively coarse grid resolution is the limited flexibility to capture multiscale features of complex bathymetry [43]. As a result, the model likely did not capture all the details of the complex bathymetry of Gnarabup Beach. Furthermore, phase-resolving models like SWASH require fine resolutions to capture higher frequency waves, e.g., [1]. However, the use of finer grid sizes in this study, which greatly increases computational cost, did not significantly improve the results, as shown in Appendix A. This may reflect that to gain the benefits of higher resolution computations requires having sufficiently high-resolution bathymetry at these fine-scales, which was not available in this study. Presently SWASH also only allows uniform wave conditions (based on single directional spectrum) to be imposed on the offshore boundary. At the study site the AWACs offshore of each reef demonstrated slightly different wave conditions which could not be accounted for in the model (which was forced uniformly with the conditions observed at site C0).
When comparing the IG waves, the model exhibited stronger IG energy dissipation and weaker nearshore currents in shallower area (channels and lagoon) as the friction increased ( Figure 6). The impact of various c f values to IG waves and mean currents has been identified in a number of reef hydrodynamic studies, e.g., [18,44,45]. In the present study, the best model performance for IG waves occurred when a uniform c f of 0.05 was applied across the entire domain ( Figure 6). While SS waves, IG waves and setup were predicted with moderate to good skill, predictions of the mean wave-driven current were often less accurate. To investigate if predictions of the currents could be improved by using a spatially-uniform c f , we considered a simulation with a spatially variable c f to investigate its impact on the skill of the modeled currents (Figure 12a). The c f values of different locations within the domain were determined from visual estimation of aerial images. Sand was assigned c f = 0.002 based on [46], 0.02 for macro-algae vegetation, e.g., [47] and 0.05 for reef, e.g., [48]. We also identified some areas dominated by a mix of sand and vegetation, in which we defined the c f value of 0.01. The offshore forcing for this simulation was the moderate wave scenario with H s = 3.4 m with a T p = 13.4 s at AWAC C0. We compared the result of bulk waves at SS-IG bands, and mean currents (magnitudes and directions) at the 22 observed points, between the spatial c f and the uniform c f scenarios. At the offshore sites, there were effectively no differences between the compared variables of both scenarios, that also suggested the bottom friction has no impact on the waves and currents at deep water. However, in shallower depths (i.e., channels and lagoon), the spatially-varying c f simulations exhibit higher magnitudes of all modeled waves heights and currents than that of the uniform c f (Figure 12b-e). Overall, using the spatially variable c f resulted in worse model performance overall than using a spatially-uniform c f (Figure 12 spectrum) to be imposed on the offshore boundary. At the study site the AWACs offshore of each reef demonstrated slightly different wave conditions which could not be accounted for in the model (which was forced uniformly with the conditions observed at site C0). When comparing the IG waves, the model exhibited stronger IG energy dissipation and weaker nearshore currents in shallower area (channels and lagoon) as the friction increased ( Figure 6). The impact of various cf values to IG waves and mean currents has been identified in a number of reef hydrodynamic studies, e.g., [18,44,45]. In the present study, the best model performance for IG waves occurred when a uniform cf of 0.05 was applied across the entire domain ( Figure 6). While SS waves, IG waves and setup were predicted with moderate to good skill, predictions of the mean wave-driven current were often less accurate. To investigate if predictions of the currents could be improved by using a spatially-uniform cf, we considered a simulation with a spatially variable cf to investigate its impact on the skill of the modeled currents (Figure 12a). The cf values of different locations within the domain were determined from visual estimation of aerial images. Sand was assigned cf = 0.002 based on [46], 0.02 for macro-algae vegetation, e.g., [47] and 0.05 for reef, e.g., [48]. We also identified some areas dominated by a mix of sand and vegetation, in which we defined the cf value of 0.01. The offshore forcing for this simulation was the moderate wave scenario with Hs = 3.4 m with a Tp = 13.4 s at AWAC C0. We compared the result of bulk waves at SS-IG bands, and mean currents (magnitudes and directions) at the 22 observed points, between the spatial cf and the uniform cf scenarios. At the offshore sites, there were effectively no differences between the compared variables of both scenarios, that also suggested the bottom friction has no impact on the waves and currents at deep water. However, in shallower depths (i.e., channels and lagoon), the spatially-varying cf simulations exhibit higher magnitudes of all modeled waves heights and currents than that of the uniform cf (Figure 12b-e). Overall, using the spatially variable cf resulted in worse model performance overall than using a spatially-uniform cf ( Figure 12). While it is possible a different set of combinations of cf values may have resulted in better performance, the computational effort in obtaining the best possible combination of cf values would likely outweigh the incremental improvement in model performance.

Conclusions
This study evaluated the performance of the non-hydrostatic phase-resolving wave flow model SWASH that was applied to a reef fronted pocket beach with complex bathymetry. The SWASH model was applied to the study site and was used to assess the performance against the observed waves, setup, and currents at 22 instrument sites that spanned the offshore, channels, and lagoon over 4 days with persistent swell ranging in height from 2.4 to 5.0 m. Despite the complex bathymetry, the model generally reproduced the observed waves (both at SS and IG bands), setup and currents throughout the domain relatively well. The largest errors were associated with the mean currents. A spatially uniform friction coefficient, c f = 0.05, was found to best reproduce the observations, which is within the range of value found in many other reef studies. An attempt to use spatially variable friction coefficients representing the different bottom types found at the site (sand, aquatic vegetation, reef) did not improve results. In general, the model performed best in areas where the most accurate bathymetry was available (closest to survey lines), indicating that in complex nearshore reef environments such as at Gnarabup Beach, very high resolution bathymetry (order 1-10 m survey resolution) is necessary to properly resolve the nearshore hydrodynamics, independent of physical processes incorporated within the model.
The headlands and rocky reefs were found to have a major influence on the nearshore hydrodynamics of Gnarabup Beach. Even though the natural structures only partially protect the beach forming semi-protected lagoon, they effectively dissipate most of the offshore SS wave energy. Near the shoreline, and in particular the area behind the reefs, the wave motions were dominated by waves in the IG band that are less influenced by the variable water depth. Moreover, wave setup within the lagoon was strongly dependent on the offshore waves. Wave energy dissipation also generated strong onshore directed currents that flowed onshore over the reefs and circulated back offshore through the channels.
SWASH output with 2 m and that of 3 m was not significant. i.e., less than 1 cm and 1 cm/s for Hs and cross-shore currents respectively, yet the error considerably increased when we used grid size ≥4 m. Hence, we considered cross-shore grid resolution of 2 or 3 m to be suitable for the study.
The optimal alongshore grid size was determined by comparing modeled H s and velocity components in 2D to a reference grid of 1 × 1 m. For this purpose, we carried out 7 simulations with similar model setting to that of the 1D model in different sizes (cross-shore size of 1-3 m and longshore size of 2-4 m, see Figure A1d) and compare the result among them. The result showed that for all the 2D grids, the errors between the tested grids and the reference of 1 × 1m for Hs and current magnitude were not higher than~0.05 m and~0.04 m/s respectively ( Figure A1b,d). We also identified grid sizes with a larger ratio between the cross-and alongshore dimension (i.e., 1 × 3 and 2 × 4 m) had larger errors. Despite relatively small differences among the tested 2D grid resolutions, there was a large difference among the simulation time needed to finish the simulations as shown in Figure A2.
The simulation with a grid size of 3 × 4 m had the lightest computation need of 137.28 core-hours for a 90 min simulation (i.e., 1 h simulation and 30 min spin-up), 3.5 times faster than the 2 × 3 m grid which we considered to be the best result. Considering the small difference of output error (with the desired accuracy), and the computational time among tested grid resolutions, we concluded the size of 3 m (cross-shore) × 4 m (alongshore) would be the optimum grid for the 2D simulation and used in this study. models is from the southern reef (red line in Figure 1a). Each simulation was run for 30 min (excluding the spin-up time of 15 min) with two equidistant vertical layers and a single wave condition of Hs of 3.14 m and Tp of 13.4 s with JONSWAP spectrum that corresponded to the average wave height during the observation period. The RMSE of modeled significant wave height ( Figure  A1a) and cross-shore velocity ( Figure A1c) parameters between 1m grid size and coarser grids showed the difference between SWASH output with 2 m and that of 3 m was not significant. i.e., less than 1 cm and 1 cm/s for Hs and cross-shore currents respectively, yet the error considerably increased when we used grid size ≥4 m. Hence, we considered cross-shore grid resolution of 2 or 3 m to be suitable for the study.
The optimal alongshore grid size was determined by comparing modeled Hs and velocity components in 2D to a reference grid of 1 × 1 m. For this purpose, we carried out 7 simulations with similar model setting to that of the 1D model in different sizes (cross-shore size of 1-3 m and longshore size of 2-4 m, see Figure A1d) and compare the result among them. The result showed that for all the 2D grids, the errors between the tested grids and the reference of 1 × 1m for Hs and current magnitude were not higher than ~0.05 m and ~0.04 m/s respectively ( Figure A1b,d). We also identified grid sizes with a larger ratio between the cross-and alongshore dimension (i.e., 1 × 3 and 2 × 4 m) had larger errors. Despite relatively small differences among the tested 2D grid resolutions, there was a large difference among the simulation time needed to finish the simulations as shown in Figure A2.
The simulation with a grid size of 3 × 4 m had the lightest computation need of 137.28 core-hours for a 90 min simulation (i.e., 1 h simulation and 30 min spin-up), 3.5 times faster than the 2 × 3 m grid which we considered to be the best result. Considering the small difference of output error (with the desired accuracy), and the computational time among tested grid resolutions, we concluded the size of 3 m (cross-shore) × 4 m (alongshore) would be the optimum grid for the 2D simulation and used in this study.