A Numerical Study of Wave-Driven Mean Flows and Setup Dynamics at a Coral Reef-Lagoon System

Two-dimensional mean wave-driven flow and setup dynamics were investigated at a reef-lagoon system at Ningaloo Reef, Western Australia, using the numerical wave-flow model, SWASH. Phase-resolved numerical simulations of the wave and flow fields, validated with highly detailed field observations (including >10 sensors through the energetic surf zone), were used to quantify the main mechanisms that govern the mean momentum balances and resulting mean current and setup patterns, with particular attention to the role of nonlinear wave shapes. Momentum balances from the phase-resolved model indicated that onshore flows near the reef crest were primarily driven by the wave force (dominated by radiation stress gradients) due to intense breaking, whereas the flow over the reef flat and inside the lagoon and channels was primarily driven by a pressure gradient. Wave setup inside the lagoon was primarily controlled by the wave force and bottom stress. The bottom stress reduced the setup on the reef flat and inside the lagoon. Excluding the bottom stress contribution in the setup balance resulted in an over prediction of the wave-setup inside the lagoon by up to 200–370%. The bottom stress was found to be caused by the combined presence of onshore directed wave-driven currents and (nonlinear) waves. Exclusion of the bottom stress contribution from nonlinear wave shapes led to an over prediction of the setup inside the lagoon by approximately 20–40%. The inclusion of the nonlinear wave shape contribution to the bottom stress term was found to be particularly relevant in reef regions that experience a net onshore mass flux over the reef crest. measurements of waves, currents and water levels with simulations from an advanced computer model to understand the physical mechanisms that determine the current patterns and water level variations at a coral reef-lagoon system in Western Australia. Friction generated by the water moving over the rough reef structures was found to reduce the mean water levels inside the lagoon. This friction was explained by the combined presence of both waves and mean currents. Furthermore, near the reef crest, the waves peak and pitch forward before they break, and this nonlinear wave shape was found to enhance the friction from the bottom. This contribution from nonlinear wave shapes however, is generally not accurately described in larger-scale computer models and should be included in such computer models to provide more reliable simulations of the water motion in coral reefs.

variations) and small scales (bottom roughness) (e.g., Lowe et al., 2010;Monismith, 2014;Symonds et al., 1995). For reef morphologies classified as barrier reefs with no coastal landmass in close proximity, strong onshore mean flows occur over the reef crest with negligible to small setup inside the lagoon (e.g., Hearn, 1999;Lowe, 2005;Sous et al., 2020). At the other morphological limit, alongshore-uniform shore-attached fringing reefs show greater similarities to open sandy beaches and may experience minimal cross-shore mass fluxes over the reef crest with maximum setup near the shoreline (e.g., Becker et al., 2014;Buckley et al., 2015;Vetter et al., 2010).
The dynamics governing mean flows and wave setup at reef sites have often been studied based on the mean (wave-averaged) momentum balance (e.g., Mei et al., 2005;Svendsen, 2005) estimated from field measurements (Arzeno et al., 2018;Vetter et al., 2010), laboratory experiments (Buckley et al., 2015(Buckley et al., , 2016Yao et al., 2017Yao et al., , 2018, or numerical models (e.g., Green et al., 2018;Lowe et al., 2009a;Sous et al., 2020;Taebi et al., 2012). Similar to the classical momentum balance at sandy beaches (e.g., Lentz & Raubenheimer, 1999;Longuet-Higgins & Stewart, 1962), intense wave breaking near the reef crest results in a wave radiation stress gradient (or wave force) that is (at least in part) balanced by a setup gradient, resulting in increased mean water levels on the reef.
Reefs are typically characterized by having large bottom roughness due to live and relic coral features and as a result, bottom stresses (drag forces) caused by the reef substrate can significantly alter the mean momentum balance (Dean & Bender, 2006;Lowe et al., 2009b). This bottom stress arises due to the presence of waves and/or currents. Nonlinear wave shapes can potentially contribute to the mean bottom stress, especially in the surf zone where waves can be highly skewed and asymmetric (Buckley et al., 2016). Due to the relatively large roughness of reefs, the bottom stress has been found to have an appreciable effect on both the mean flow and setup dynamics in reef environments (Monismith, 2014;Symonds et al., 1995).
The impact of mean bottom stresses, however, highly depends on the reef geometry and the associated flow. Fringing reefs with large channels typically develop circulation cells with onshore flow over the reef crest and offshore flow through the channel (e.g., Lowe et al., 2015Lowe et al., , 2009bMonismith, 2014), analogous to rip currents on barred beaches (e.g., Dalrymple et al., 2011;MacMahan et al., 2006). In the case of onshore flows across the reef, the mean bottom stress generated by these flows and/or waves reduces the setup inside the lagoon and near the shoreline (Hearn, 1999;Sous et al., 2020). In contrast, for alongshore uniform fringing reefs lacking channels with zero net mass flux, a wave-driven onshore mass flux near the surface is balanced by offshore directed near-bed flows (e.g., Buckley et al., 2015Buckley et al., , 2016, similar to the undertow profile at sandy beaches (e.g., Faria et al., 2000). In this case, the bottom stress increases wave setup and the elevated water level is nearly constant landward of the surf zone, resulting in a maximum setup near the shoreline (Apotsos et al., 2007;Buckley et al., 2016).
Most detailed studies of momentum balances on reefs have been based on laboratory experiments with (simplified) reef morphologies. Albeit these studies were typically conducted with high spatial resolution, they typically focused on cross-shore dynamics and did not capture the complexity of the two-dimensional hydrodynamics that occur at natural reefs. Such two-dimensional dynamics have been measured in the field, but field experiments were often limited by sparse measurements with minimal to no sensors in the surf zone region (where critical wave transformations occur). This limits our understanding of the momentum balance at natural reef sites and complicates a rigorous assessment of the various numerical modeling approaches that are available to simulate reef hydrodynamics.
To date, most field-scale numerical studies of reef hydrodynamics have used circulation models coupled to phase-averaged wave models (e.g., Hoeke et al., 2013;Lowe et al., 2009a;Van Dongeren et al., 2013). This approach can describe both barotropic and baroclinic flows but does not intrinsically capture the nonlinearity of the wavefield and instead relies on parametrizations of their effects (e.g., Ruessink et al., 2012). In contrast, phase-resolving models, such as Boussinesq-type or nonhydrostatic models, can intrinsically capture the nonlinearity of the wavefield. Due to their increased computational requirements however, they often focused on a single cross-shore reef transect (Beetham et al., 2016;Nwogu & Demirbilek, 2010;Sous et al., 2019) and have been seldom applied at the scale of a complete reef-lagoon region (Roeber & Bricker, 2015). As such, the contribution of nonlinear wave shapes and energy transfers to the two-dimensional setup and mean flow dynamics in realistic reef regions remains largely unknown.
RIJNSDORP ET AL. In this work, we combine a highly spatially resolved field data set with a phase-resolving wave-flow model to study the mean flow and setup balance at Ningaloo Reef, located in the northwest of Western Australia. The reef system is characterized by a steep fore-reef slope, which is backed by a lagoon that is connected to the ocean through channels. In this region, flows are predominantly wave-induced with minor contributions from wind and tides, and negligible buoyancy effects (Taebi et al., 2011). Measurements were obtained at an unprecedented scale with >15 sensors across the surf zone at the reef crest (where the strongest gradients in the waves, flow, and setup occur), allowing for a rigorous assessment of the capability of the model in reproducing the local hydrodynamics. The hydrodynamics at the site (spanning approximately  4 6 km) were modeled using the nonhydrostatic wave-flow model SWASH (Zijlema et al., 2011). Model predictions were used to quantify the main drivers of the mean flow and setup dynamics through the mean momentum balance with particular attention to the role of nonlinear wave shapes.
In Section 2, we present the methodology used in this paper, including a description of the field experiment, the numerical model and its setup, and the data analysis methods. Subsequently, we compare model predictions and field observations for two consecutive swell events occurring over a 2-day period (Section 3). This comparison showed that the model captured the wave transformation and the prevailing flow patterns in this complex reef system, allowing for a more detailed analysis of the modeled mean momentum balance to study the mean flow and setup (Section 4). First, we study the cross-shore momentum balance to understand the main drivers of the mean flow and setup dynamics (Section 4.1). This is followed by a decomposition of the bottom stress term into its mean flow and (linear) wave contributions in order to understand how these contributions alter the setup inside the lagoon (Section 4.2). Detailed cross-shore simulations over the central reef transect were conducted to support the modeling at the scale of the reef-lagoon system (Section 5). The main findings of this work are summarized in Section 6.

Field Observations
The field study was conducted between May and June 2016 at a section of the Ningaloo Reef, Western Australia (21°52′9.66″S, 113°59′26.95″E). The site is characterized by a steep 1/20 fore reef slope, a relatively narrow (100s of m) approximately 1 m deep reef flat, and a 2 km wide shallow lagoon (<5 m deep). Inside the lagoon the seabed is predominately sandy, whereas the seabed around the reef crest is much rougher due to relic and living coral . The site is micro-tidal with an average tidal range of 0.8 m ( Figure 1e). During the experiment, waves predominantly arrived from the southwest, with typical wave heights 0 m H of 0.5-2 m and peak wave periods T p of 10-18 s ( Figure 1d). During extreme weather conditions, local winds can contribute significantly to the hydrodynamics in the lagoon (e.g., Cuttler et al., 2018;Drost et al., 2019). During the study period, the wind was typically weak to moderately strong, varying around an average of approximately 5 m/s (Figure 1c), and its contribution to local wave generation and setup in the lagoon was estimated to be small (see Section 2.3, for further details). The tidal component of the currents was found to be typically small in this region (Taebi et al., 2011), which was confirmed by the strong correlation of the currents with offshore wave height (see Section 3).
A combination of velocity and pressure sensors was deployed at the field site (see Figure 1a) during the approximately 6 weeks experiment (see Table 1, for an overview and settings of the instruments). In particular, a dense cross-shore transect of pressure sensors was placed over the reef crest to measure the wave and setup dynamics throughout the surf zone. Pressure sensors were deployed on hard substrate to prevent settling, and their measurements were converted to water level signals using linear wave theory. At all pressure sensors, we computed the setup as the change in the mean water level relative to the deepest pressure sensor (located in approximately 18 m water depth) following the approach of Lowe et al. (2009b), in which   0 ΔP P P is the difference between the (hourly averaged) water depth (assuming hydrostatic pressure) at the respective sensor and the reference sensor (indicated by subscript 0), and h is the still water depth (excluding any tidal, wind and wave contributions). The still water depth was estimated at instances RIJNSDORP ET AL.

10.1029/2020JC016811
3 of 22 when negligible setup was expected (i.e., in the case of relatively small waves with short periods, low wind speeds, and near high tide). We subsequently used the still-water depth estimate closest to the respective burst to compute the setup estimate from Equation 1. The measurements from the various devices were analyzed in bursts with a duration of 1 h.
Simulating the entire 6 weeks was not feasible due to the computational requirements of the phase-resolving wave-flow model (simulating an individual sea-state at this reef-lagoon system for a duration of ∼1 h takes up to 24 h on ∼100 cores at the Pawsey Supercomputing Center where the simulations were conducted). We therefore modeled a subset of the conditions measured during the field campaign, which captured two swell events that occurred between the May 30 and the June 3 (Figures 1b and 1c). A total of 35 individual wave conditions were modeled (once every approximately 3 h, resulting in a simulation near high, mid, and low tide during each tidal cycle).

Numerical Model, Governing Equations
We used the three-dimensional nonhydrostatic wave-flow model SWASH (Zijlema et al., 2011) to simulate the wave conditions, which solves the Reynolds-Averaged Navier-Stokes equations. The model considers a three-dimensional fluid domain that is bounded in the vertical by the free-surface x y t and bottom     , z d x y , in which t is time and x,y,z are the Cartesian coordinates (with  0 z at the still water level). In this framework, the governing equations read, using the Einstein summation convention for the horizontal coordinates, where u i is the horizontal component of the velocity in the i direction, w is the vertical velocity component, ζ is the free-surface, p is the total pressure (hydrostatic plus nonhydrostatic pressure), g is the gravitational constant, and τ represents the turbulent stresses.
A quadratic friction law was used to model the stress induced by the sandy bottom and the reef structures at the seabed, in which c f is a friction coefficient, and U and V represent the depth-averaged velocity in the x and y direction, respectively. The bottom stress components are given by    , where in SWASH the angle θ is taken from the horizontal velocities in the bottom layer (u b and v b , respectively) As the magnitude of the bottom stress depends on depth-averaged velocities, it does not take into account the vertical structure of the mean flow, the frequency-dependent vertical attenuation of the orbital velocities of the waves, nor the flow within large bottom roughness features (a canopy layer) of the reef. This is a limitation of the lagoon-scale modeling presented in this study (Sections 3 and 4).
To confirm the findings of the lagoon-scale modeling, we therefore conducted additional modeling with a more sophisticated description of the rough reef that does capture this frequency-dependent behavior (Sections 2.4 and 5.1).
For the turbulent terms, the model uses separate eddy viscosities to account for horizontal and vertical mixing (for further details see Rijnsdorp et al., 2017). In this work, the horizontal eddy viscosity was estimated based on a Smagorinsky-type formulation (with default parameters), in combination with a small and constant vertical eddy viscosity to enhance vertical mixing (   5 1 10 ). To capture the onset of wave breaking, we used the Hydrostatic Front Approximation with default parameters (Smit et al., 2013).

Numerical Model Setup-Lagoon Scale Modeling
The bathymetric grid was constructed from gridded 0.25 m resolution LiDAR data (collected by Airborne Research Australia using a Riegel Bathymetric Lidar), single beam DGPS-corrected bathymetric surveys, and gridded 50 m resolution multibeam data (Wilson et al., 2012). LiDAR data were used for 0-20 m water depths, and multibeam data for water depths greater than 20 m. Single beam DGPS-corrected bathymetric RIJNSDORP ET AL.
10.1029/2020JC016811 5 of 22 surveys were conducted in order to fill gaps between the coverage of the LiDAR data and the multibeam data. An artificial beach profile was used upward from the mean water level (with the same 1/20 slope along the whole domain), which preserved the shape of the actual shoreline (including the salient in the center of the lagoon). This was done to replicate the typical beach profile, while avoiding issues related to short-term spatially varying beach profile changes. The maximum still water depth (excluding tidal level) in the model domain was set at 18 m, which is equal to the depth at the offshore AWAC used to force the numerical wavemaker.
The alongshore variability of the reef region required a careful selection of the numerical domain to capture the dominant wave and mean flow dynamics across the study area. To reduce the required domain size and keep the computational requirements feasible, we used cyclic boundary conditions to reduce boundary effects at the lateral (northern/southern) boundaries. The cyclic boundary implies that waves and flow exiting through one boundary reenter the domain through the opposite boundary. The selected model domain spanned a region that included the main reef of interest, and the main channels at the northern and southern reef edges ( Figure 1a). To prevent discontinuities in the bathymetry at the lateral boundaries (due to the use of cyclic boundary conditions), the measured bathymetry was extended with a 200 m wide linear transition between the two outermost cross-shore transects. Although this domain does not include the natural bathymetry toward the north and south of the reef-lagoon system, we found that extending the domain to capture both the exposed coast to the north and the lagoon to the south did not significantly affect the wave and mean flow fields in the region of interest (not shown).
The final numerical domain spanned  4,590 6,560 m, including a 100-200 m wide constant depth region near the wavemaker. For each simulation, the still water level was offset with a constant value according to the measured tidal level. We did not attempt to capture the tidal currents within the domain as it is currently not possible to incorporate tidal currents in SWASH. To confirm that tidal velocities are small at this field site, we conducted additional simulations (tidal flow only) using the Delft3D-Flow model setup of Cuttler et al. (2018). Results of these model simulations showed that tidal velocities were <0.1 m/s (including near the reef crest and at the channels) during the subset of conditions that were considered in this work (not shown). This confirms that tidal velocities were typically smaller than the wave-induced velocities (which varied between 0.1 and 0.5 m/s, see e.g., Figure 5). Furthermore, the contribution from locally wind-generated waves to the sea-swell waves within the reef was estimated to be negligible small. Using the wave growth relation of Breugem and Holthuijsen (2007), with a fetch of 2 km and depth of 5 m (based on the geometry of the reef-lagoon system) and a wind speed of U 10 = 3-7 m/s (based on the measured local wind conditions), we estimated that the locally generated wind waves in the lagoon had a significant wave height of The locally generated wind-wave energy thus resides at frequencies that fall outside the sea-swell frequency band (defined in Section 2.5). The wind contribution to the setup in the lagoon was also estimated to be small. Using the nonlinear shallow water equations with a wind stress at the sea-surface, the wind-setup at the shoreline for an onshore directed U 10 = 3-7 m/s was estimated to be     Ο 0.01 m. As such, we did not attempt to include their contributions in the SWASH simulation.
At the offshore boundary of the domain, a directional sea-state was generated based on measured two-dimensional wave-spectra from the offshore AWAC. Wavefields were generated using linear wave theory and a second order correction to include bound infragravity waves (Rijnsdorp et al., 2014(Rijnsdorp et al., , 2015. At the shoreline, a moving-shoreline boundary condition was used to capture the wetting and drying of the beach (Stelling & Duinmeijer, 2003).
The model domain was discretized using a regular grid with a horizontal resolution of  Δ 3 x m and  Δ 5 y m, which was found to be sufficiently fine based on a grid convergence study (Appendix A1). Two vertical layers were used to capture the dispersive properties of the sea-swell waves. To ensure model stability, the time step was changed dynamically during the simulation to ensure that the CFL number did not exceed 0.8 (with an initial time step of  Δ 0.1 t s). For each simulation, the model was run for 90 min. This includes 30 min spin-up time to ensure that the total kinetic energy, potential energy, and entropy reached a quasi-steady equilibrium (following Feddersen et al., 2011). Water level and layer-averaged velocity signals were outputted for the last 60 min of the simulation at 2 Hz (similar to the measurements) throughout the domain.
10.1029/2020JC016811 6 of 22 The bottom friction coefficient was set to vary inside the domain, with  0.005 f c for the sandy regions and with enhanced frictional values near the reef crest to mimic the frictional dissipation induced by the reef structures. An increased friction coefficient was imposed over the fore-reef, reef crest and reef flat according to a calibration study for three representative sea-states out of the total of 35 wave conditions that were considered (spanning relatively large, moderate, and small wave heights). Overall best agreement between model predictions and observations of the setup and bulk wave height (sea-swell and infragravity) was found for  0.02 f c (Appendix A2). This friction coefficient falls within the range previously reported for coral reefs (Lowe et al., 2009a;Rosman & Hench, 2011).

Numerical Model Setup-Cross-Shore Modeling Over the Central Reef Transect
Due to the computational limitations associated with the spatial scale of the reef-lagoon system, the lagoon-scale modeling was restricted to a coarse vertical resolution (two vertical layers were used in the simulations) in combination with a simple quadratic friction law to describe the bottom stress. As a result of the coarse vertical resolution, the lagoon-scale modeling does not account for the vertical structure of the mean flow, the (frequency dependent) vertical attenuation of orbital velocities, and the flow inside the canopy-like reef structures. To confirm the findings of the lagoon-scale modeling, we performed an additional set of SWASH simulations to assess the influence of a more detailed description of the vertical flow structure and the interactions of the flow with the rough reef. The rough reef was described using a canopy (or vegetation) model (Suzuki et al., 2019) that takes into account the vertical structure of the mean flow and orbital velocities, and resolves the flow inside a (simplified) canopy. This approach provides a more accurate description of the rough reef compared to the quadratic friction flaw.
To account for the draining of the lagoon through the channels that (partially) compensated for the onshore wave-induced mass flux over the reef, we implemented a recirculating pump system within SWASH (conceptually similar to the laboratory experiment of Yao et al., 2018). The discharge through the pump was driven by the pressure head difference between the lagoon and the offshore region (seaward of the reef crest) through the Bernoulli principle. For a unit width, the resulting discharge was computed as, in which η L and η o are the instantaneous setup (computed as a moving average over a 5 min duration) inside the lagoon and offshore of the reef crest, respectively. Including this discharge resulted in a near complete compensation of the onshore wave-induced mass flux and a comparable mean water level inside and outside of the lagoon, which is representative of a barrier type reef system (e.g., Sous et al., 2020). However, frictional dissipation may limit the effectiveness of the outflow through the channel, and many reef-lagoon systems are characterized by a water level difference between inside and outside of the lagoon (including the reef-lagoon system considered in this work). We therefore carried out different simulations in which a fraction of the full discharge was considered (αQ pump , with α ranging between 0 and 1). This allowed us to cover the natural variability of reef systems, ranging from fringing reefs with no back-lagoon and channels (  0) to barrier reefs (  1).
All previously considered 35 wave conditions were simulated with the cross-shore two-dimensional vertical (2DV) SWASH model for    0 1 with 0.25 increments, resulting in a total of 175 simulations. Simulations were performed with 20 layers, and with the same horizontal resolution and initial time step as the two-layer simulations of the whole reef-lagoon system. The total duration of the simulations was 90 min, including 30 min spin-up. A single 2DV simulation took   Ο 1h to complete on an 8-core desktop machine. The  k  turbulence model with default parameters was used to capture the vertical turbulent terms. At the bottom boundary, the law of the wall was used to determine the near-bed velocities, and to prescribe  The rough reef structures were included using a canopy flow model (Suzuki et al., 2019), with which we idealize the rough reef as a canopy composed of vertical cylinders that interact with the flow (e.g., Lowe et al., 2005). With this approach, the force on the flow exerted by vertical cylinders (f c ) was included in the momentum equations based on the Morison equation (Morison et al., 1950), Apart from the inclusion of an inertia term (the second term in the right-hand side of Equation 8), the main difference compared to the quadratic friction law, Equation 6, is that the canopy force f c depends on the local flow velocities inside the canopy, whereas the quadratic friction law depends on the depth-averaged flow velocities. This dependence on the local instead of the depth-averaged velocities can have important implications, as this approach captures the frequency-dependent dissipation of wave energy, allowing for a more accurate representation of dissipation by canopies (e.g., Jacobsen et al., 2019;Jadhav et al., 2013;Lowe et al., 2007). The canopy parameters in this equation, their physical meaning, and their value as used in this study are displayed in Table 2. These parameters were chosen such that the combined value of all parameters was equal to the constant friction coefficient of the quadratic drag law used in the two-layer simulations (i.e.,  0.02 f c ).

Data Analysis
To facilitate a model-data comparison, predicted and measured water level and velocity signals were processed in the same manner. Wave spectra (E f ) were computed from the water level signals with 30 degrees of freedom. Spectra were divided into a sea-swell (SS) band 2 p p f f f based on the peak period f p of the incident wave spectrum (Janssen et al., 2003;Pomeroy et al., 2012). Bulk significant wave heights were computed for both bands by integrating the spectra over their respective frequencies e.g., H E df Measured velocity signals were rotated into the local coordinate system that was used in the model (consistent with the axis in Figure 1a) and vertically averaged in the case of the ADCP instruments. Depth-averaged velocities were output from the model. The resulting measured and predicted velocities were subsequently averaged in time (over 60 min) for each sea-state. We thus compared the predicted depth-averaged velocities with both point measurements (at the ADV instruments) and vertically averaged velocity profiles (at the ADCP instruments).
To quantify the model performance, we computed the relative bias, RB (e.g., Van der Westhuysen, 2010) and Murphy's skill score, MS (Murphy, 1988) at each instrument and for all sea-states considered, where Q is a predicted (subscript P) or measured (M) quantity, and N is the total number of measured/ predicted values. The relative bias indicates whether the model underpredicts (negative RB) or overpredicts (positive RB) the observations. The skill score is equal to one for perfect agreement and a score below zero indicates that the predictive ability of the model is worse than assuming the mean of the measurements.

Mean Momentum Balance
To understand the main driving forces of the (wave-driven) mean flow and setup, we investigated the wave-averaged and depth-averaged momentum balance terms computed from the model (e.g., Mei et al., 2005). The different momentum terms were computed during the simulations using the implementation of da Silva et al. (2020), which greatly reduced data requirements and ensured the computed momentum terms are consistent with the numerical discretization of the model equations (resulting in a residual that is in the order of machine precision). In a Cartesian framework, the modeled mean momentum balance reads, in which the brackets  denote time averaging,      h d is the mean water depth with η equal to the mean water level or setup (   ), U L,i is the depth-averaged mass flux velocity (or Lagrangian flow velocity, e.g., Rogers et al., 2013) in which   z d q is the nonhydrostatic pressure at the bottom, and S ij represents the radiation stresses. The radiation stresses S ij are given by, in which δ ij is the delta Dirac function, and  i u are the oscillatory velocities in the i direction (      , i L i u u U ). As the model solves the RANS equations, the oscillatory velocities do not include any turbulent fluctuations, and S ij can thus be identified as the wave-radiation stress.

Bottom Stress Decomposition
The mean bottom stress term in Equation 9 integrates the total contributions of mean currents (nonlinear) wave velocities, and interactions between the two. To differentiate between these contributions, we separated the bottom stress from Equation 6 into A mean flow, wave, and a wave-current interaction contribution (e.g., Rooijen et al., 2020). For this purpose, we decomposed the velocity signals in a mean and oscillatory component (e.g.,      U U U). For example, substituting this relation into the x-component of the bottom stress parametrization, and averaging in time yields, in which the first term   The oscillatory velocities include the potential contribution from nonlinear wave shapes, as the signals can have nonzero skewness and asymmetry (especially in the surf zone region near the reef crest). To evaluate the contribution from the nonlinear wave shapes to the bottom stress, we also estimated the mean bottom stress (Equation 14) for a linearized wavefield with zero asymmetry and skewness. For this purpose, we transformed each velocity signal used in the bottom stress parametrization ( , , b U V u , and v b ) to the frequency domain (using the Fourier transform), and added a random phase to all Fourier components. The same random phases were added for each set of Fourier components throughout the numerical domain, which ensured that the signals maintained their phase dependency in space (e.g., maintaining the wave directionality). The addition of random phases to the frequency components removed any phase locking between primary waves and their higher order harmonics. Transforming the resulting signals back to the time-domain using the inverse Fourier transform resulted in velocity signals with zero skewness and asymmetry, which can be interpreted as a linear representation of the original velocity signals. Computing the bottom stress based on the resulting linearized wave signals (including the mean flow contribution) yielded an estimate of the bottom stress due to the combination of mean flow and linear waves. In the following, we refer to the total mean bottom stress from linear velocity signals as, is the mean stress due to the interaction between the linear waves and the mean flow.

Model-Data Comparison
The model captured the typical bulk wave heights at the sea-swell frequencies (H SS , Figure 2a) with overall good skill and small bias (  0.97 MS and  0.05 RB , see Figure 3a). This includes the increase of H SS as sea-swell waves shoal over the fore-reef slope toward the reef crest, and the depth-limited reduction of H SS (including tidal modulation) at sensors near the reef crest due to intense wave-breaking on the shallow reef (Figures 2a and 4a). Predicted and observed infragravity wave heights H IG were generally in agreement (  0.72 MS and  0.02 RB , Figure 3b), with H IG peaking near the reef crest and gradually reducing toward the lagoon, where H IG was typically slightly larger than at the offshore sensor (Figures 2b and 4a). The model also reproduced the cross-shore variation of wave setup η (  0.89 MS , Figure 3c). This includes the set-down just prior to wave breaking, an increased η over the reef crest (albeit occasional under-predictions, as illustrated by the negative bias,  0.22 RB ), and smaller η inside the lagoon than at the reef crest (Figures 2c and 4b).
Observed mean flow velocities, in the local coordinate system oriented with the reef, were indicative of a general circulation pattern with strong shoreward flow over the reef crest and diverging flows directed toward the channels inside the lagoon (Figures 4c and 5a). Flows were typically largest over the reef crest and in the channels and were smaller inside the lagoon. Modeled and observed net velocity vectors and velocity variance ellipses (computed over the range of modeled wave conditions) show that the model captured the dominant circulation patterns and its variability in the reef system (Figures 5a and 5b). The velocity sensor near the tip of the salient (at  2 x km,  0 y km) forms a notable exception. Although the velocity magnitude and variance were comparable, the predicted velocity was directed toward the northern channel whereas the measured velocity was directed toward the southern channel. This discrepancy between the predicted and measured flow direction is likely related to the sensor being located in the region where the onshore flow diverges toward northern/southern direction. In this region, the flow direction changes over relatively small spatial scales, and as a result, measurements of the flow direction are sensitive to small deviations in the sensor location. The generally good model-data agreement is further exemplified by scatter plots of the mean flow velocities versus the offshore wave height (Figures 5c and 5d). The good agreement shows that the model captured the increasing mean flow velocities for larger incident wave heights, although velocities are somewhat under predicted at the ADV instruments near the reef crest (e.g., S1 and Figure 5c) and over predicted at the ADCP profilers located near the edge of the reef flat (e.g., S2 and Figure 5d). At RIJNSDORP ET AL.    Figure 3d). Cross-shore velocities were typically larger than alongshore velocities, resulting in a good reproduction of the velocity vector magnitude (  0.72 MS , Figure 3f).

Cross-Shore Mean Momentum Balance
To identify the main drivers of the setup dynamics, we studied the contribution of the different terms to the modeled cross-shore mean momentum balance. In the following, we will first inspect the balance in detail for the most energetic sea-state (the first wave condition considered, see Figure 1), followed by a quantification of how the different terms contributed to the setup balance for all wave conditions considered. This allows us to identify the dominant drivers of the setup dynamics inside the reef-lagoon system.
The dominant terms that contributed to the modeled mean cross-shore momentum balance at three different reef transects during the most energetic sea-state are shown in Figures 6b-6d. Over the central reef transect (T3, Figure 6c), wave breaking near the reef crest resulted in a large wave force W f,x that was primarily balanced by a cross-shore pressure gradient ). The imbalance between these two terms provided an excess force that was primarily balanced by a combination of the bottom stress τ b,x and the advection terms    , , L j x L i j h U U , with τ b,x providing the largest contribution. The various terms typically peaked near the reef crest, and gradually reduced over the reef flat toward the lagoon. The excess forcing (net of pressure gradient and wave force) remained nonzero up to  x 500 m (approximately 500 m shoreward of the reef crest), where the contributions from the wave force and the advection terms were relatively small. In this reef flat region    ( 100 m 500 m) x , a negative pressure gradient was primarily balanced by the bottom stress (similar to the observations of Sous et al., 2020 for a barrier reef). These patterns were consistent across the length of the reef, as illustrated by the two transects (T1 and T5) located in close vicinity to the southern and northern channel (Figures 6b and 6d, respectively). Inside the surf zone, the net excess forcing from the imbalance between the wave force and pressure RIJNSDORP ET AL.  gradient was compensated by primarily the bottom stress with a smaller contribution from the advection terms. On the reef flat, the excess forcing was also primarily balanced by the bottom stress, but in contrast to the surf zone, the excess forcing was primarily due to the pressure gradient, given that the wave force was small on the reef flat.
To quantify how the different momentum terms contributed to the setup inside the lagoon, we reconstructed the setup profile by numerically integrating the mean cross-shore momentum balance Equation 9 for an unknown setup η. By including different combinations of forcing terms, we estimated the influence of the various forcing terms on the setup balance. For example, the contribution from the wave force and bottom stress can be evaluated by numerically integrating the setup based on the following balance (using a simple finite difference approach), We compare the reconstructed setup inside the lagoon η L (taken 500 m from the shoreline) based on a partial momentum balance (vertical axis) and the full momentum balance that includes all momentum terms (horizontal axis) in Figure 7. To gain insight into how the main drivers of the setup balance varied along the lagoon, we distinguished between five regions (centered around the transects plotted in Figure 6) that span the full alongshore width of the reef-lagoon system. Although the pressure gradient and wave force are approximately of equal magnitude (Figure 6), reconstructing the setup solely based on W f,x led to a significant over prediction (  2 RB , indicative of typical over predictions that exceed 200%) of the setup inside the lagoon (η L ) for all five regions (Figure 7). Inclusion of the bottom stress term greatly reduced η L and improved the accuracy of the reconstructed η L , resulting in a (near) zero RB (i.e., no consistent over or under prediction) but still some scatter. Toward the channels, the scatter in this reconstructed η L from Equation 16 typically increased, as indicated by the lower skill scores relative to the central reef transect. The inclusion of the advective terms was of increasing importance toward the channels (and in particular in region T5, RIJNSDORP ET AL.

Bottom Stress Decomposition
Thus far, the bottom stress term has incorporated all possible contributions related to the mean flow, the (nonlinear) wave velocities, and interactions between the mean flow and the waves. To differentiate between these contributions, we diagnostically evaluated the bottom stress from Equation 11 using three different velocity signals as outlined in Section 2.7. Figure 8a evaluates how the (cross-shore integrated) "reconstructed" bottom stress computed based on the three velocity signals compares to the total mean bottom stress for all wave conditions and along the whole reef-lagoon system. represents the cross-shore pressure gradient, W f,x the cross-shore wave force, 〈τ b,x 〉 the mean cross-shore bottom stress, and    , , L j x L i j h U U the advective contributions. The pressure gradient and wave force are the largest contributors, and values larger than the axis limits are not shown. The local still water depth d (measured positive down) is plotted below each panel. For reference, panel (a) shows the predicted setup throughout the numerical domain (colors with the bathymetry indicated by the contour lines), the measured setup (colored circles), and the location of the considered transects (horizontal dashed black lines). These transects also define the center of the five regions that will be considered in Figure 7, with their boundaries indicated by the full black lines. nonlinear wave shape was required to obtain perfect agreement and to fully explain the total bottom stress (not shown).
To quantitatively assess how neglecting parts of the velocity signal impacted the setup inside the lagoon, we repeat the analysis of Figure 7 including all momentum terms, but with the reconstructed bottom stress term (from all three methods) instead of the actual modeled bottom stress (which was used in the previous section). Consistent  indicates that the exclusion of the contribution from nonlinear wave shapes to the bottom stress resulted in an average over prediction of the setup inside the lagoon by approximately 40%.

Discussion
A phase-resolving nonhydrostatic wave-flow model (SWASH) was used to study the mean flow and setup balance at a section of Ningaloo Reef, Western Australia. Modeled bulk wave parameters (wave heights and setup) and mean flow velocities were compared with extensive field measurements for a time span of approximately 4 days (covering two consecutive swell events). The model reproduced the sea-swell wave height and setup in the reef system with excellent skill (0.97 and 0.89, respectively). It captured the wave breaking over the shallow reef crest, the saturated (tidally modulated) sea-swell wavefield on the reef flat and inside the lagoon, and the cross-shore setup structure with maximum setup near the reef crest and a reduced setup of up to ∼50% inside the lagoon (relative to the reef crest). The model also captured the RIJNSDORP ET AL.  Comparison between the reconstructed setup inside the lagoon (η L ) based on the partial (various terms excluded) mean cross-shore momentum balance (vertical axis) and η L reconstructed using all mean cross-shore momentum terms (horizontal axis). The five subplots represent a subregion of the reef-lagoon system, which combined span the full reef-lagoon system (as indicated in Figure 6a). The setup in the lagoon was taken 500 m from the shoreline, at all alongshore locations within the five subregions. The setup was reconstructed using a combination of mean cross-shore momentum terms, as indicated by the marker color and the legend in the bottom right (in which, W f,x represents the cross-shore wave force, 〈τ b,x 〉 the mean cross-shore bottom stress, and    , , L j x L i j h U U the advective contributions). Each panel includes a legend with the skill score (MS) and relative bias (RB) of each approach to reconstruct the setup.
patterns and variability of the mean (wave-averaged) flow, with onshore directed flows over the reef crest into the lagoon and concentrated outflow through the channels. Discrepancies between observed and modeled mean flows were typically larger than for the bulk wave parameters, with typical under predictions of the alongshore current component near the reef crest (resulting in a poor skill of 0.22). Part of these discrepancies could be caused by the absence of tidal flows in the SWASH simulations. Nonetheless, the model reproduced the dominant flow patterns and its variability and captured the onshore flow components and absolute current magnitude with reasonable skill (0.77 and 0.72, respectively).
The modeled cross-shore mean momentum balance indicated that the wave force was primarily balanced by the pressure gradient and bottom stress over the central reef section. Advective terms became increasingly important toward the channels at the northern and southern end of the reef system. On the reef flat, located shoreward of the surf zone, the flow was primarily driven by a setup gradient as the wave force was generally small. Associated with an elevated setup inside the lagoon relative to offshore, the pressure gradient also provided the main forcing of the flow inside the lagoon, driving strong seaward flows through the channels (not shown).
The mean momentum balances further highlighted that the setup inside the lagoon was primarily controlled by the wave force and the bottom stress. The bottom stress induced by the rough reef resulted in an appreciable reduction of the setup inside the lagoon. The bottom stress was primarily explained by the combined effect of waves and currents. Excluding the contribution from nonlinear wave shapes resulted in an appreciable over prediction (∼40%) of the setup inside the lagoon.

Implications of Cross-Reef Mass Flux and the Vertical Flow Structure on the Setup Profile: Additional 2DV Modeling
Due to the computational limitations associated with the spatial scale of the reef-lagoon system, the modeling was restricted to a coarse vertical resolution (two vertical layers were used in the simulations) RIJNSDORP ET AL.
10.1029/2020JC016811 16 of 22 from the nonlinear wave shapes of the mean canopy force was required to remove this bias and accurately reconstruct η L for all α values (not shown).
The 2DV modeling did not account for wave directionality, and the complexity of real reef sites (such as the alongshore reef variability and the presence of channels). Nonetheless, the introduction of the onshore mass flux through the pump system in the model allowed for the inclusion of the effect of the draining of a reef flat or lagoon (e.g., through channels) on the cross-shore reef hydrodynamics. We therefore believe that the results are representative for the hydrodynamics at a centrally located reef transect (e.g., sufficiently away from the channel) where the dominant wave motions and wave-induced mean flows are predominantly cross-shore directed (e.g., Figure 5). The results support the findings of the full-scale modeling resolving the 2D bathymetry and indicate that nonlinear wave shapes can contribute significantly to the mean bottom stress; when excluding the effect of nonlinear wave shape, setup inside the lagoon was overpredicted by up to 20%. The inclusion of this bottom stress contribution was found to be particularly important for reefs that experience a net onshore mass, as the mean bottom stress did not significantly alter the reconstructed η L for   0 (Table 3). This finding is particularly relevant for future phase-averaged modeling efforts, which require explicit parametrization of nonlinear wave shapes (e.g., Ruessink et al., 2012) to capture their contribution to the mean bottom stress (e.g., van Rooijen et al., 2016).

Conclusions
Detailed field observations and phase-resolved numerical simulations (with the non-hydrostatic wave-flow model, SWASH) were used to investigate the mean flow dynamics and setup balances over a reef-lagoon system (Ningaloo Reef, Western Australia). Strong onshore flows over the reef crest into the lagoon were shown to be driven by the wave force (dominated by gradients in radiation stress) due to breaking at the reef crest. The flow diverged inside the lagoon and was directed toward the channels, where strong outflows drained the lagoon. The cross-shore setup balance varied at different alongshore transects of the reef system. Over the central reef transect, the setup inside the lagoon was primarily controlled by the wave force and bottom stress. Toward the channels at the lateral edges of the reef system, the advection term became increasingly important, but excluding it when reconstructing the setup inside the lagoon resulted RIJNSDORP ET AL.
in increased scatter (but no bias). The bottom stress, generated by the interaction of the waves and flow with the relatively rough reef, acted to significantly reduce the setup inside the lagoon. Excluding the effect of bottom friction in the setup-balance resulted in an over prediction of the wave-setup inside the lagoon of 200-370%. The bottom stress was explained by the combined presence of waves and wave-induced currents. The contribution from the wave and wave-flow interactions were essential to fully explain the modeled bottom stress. Nonlinear wave shapes resulted in an increase of the bottom stress (especially near the reef crest where wave breaking occurs) and exclusion of their contribution resulted in an over prediction of the wave-setup inside the lagoon of up to 20-40%. Additional modeling suggested that the contribution of nonlinear wave shapes to the bottom stress is particularly relevant for reef systems that experience an onshore mass flux over the reef crest.

GRID Convergence
A grid convergence study was conducted to decide on the optimal grid resolution for the two-layer simulations at the scale of the reef-lagoon system. A single energetic sea-state was simulated with variable grid resolutions (   Δ 2 10 x m and   Δ 2.5 10 y m), and we computed the skill score of various bulk wave and flow parameters for each simulation relative to the finest grid resolution (

Friction Coefficient
The enlarged friction coefficient over the rough fore-reef and reef flat was chosen based on a calibration study. Three representative wave conditions (covering a relative energetic, moderate, and mild sea-state) were simulated with different c f values over the fore-reef and reef flat. For each simulation, we computed the skill score in reproducing the bulk wave heights (sea-swell and infragravity) and setup that were measured at the instrument sites. Best agreement for all three parameters was found for  0.02 f c (Figures 11b), which was subsequently used in all two-layer simulations.

Data Availability Statement
The processed results of the field experiment that are displayed in the figures, and the input files for the numerical model are publicly available at https://doi.org/10.26182/wt45-eb58. RIJNSDORP ET AL.  ). Panel (b) shows the skill score (relative to the measurements) for the sensors over the main reef crest for three wave conditions (with a relatively energetic, moderate and weak sea-state) as a function of different friction coefficients over the reef crest.