PhD Edward Paul Beetham
email: e.beetham@auckland.ac.nz

PhD Edward Paul Beetham

PhD Edward Paul Beetham

PhD Paul Simon Kench

PhD Stéphane Popinet

Ryan Lowe
email: ryan.lowe@uwa.edu.au

Robert Brander
email: rbrander@unsw.edu.au

Eddie Beetham

Edward Beetham

Paul S Kench

Beetham

Popinet Kench

Model Skill and Sensitivity for Simulating Wave Processes on Coral Reefs Using a Shock-Capturing Green-Naghdi Solver

Keywords: wave transformation, Green-Naghdi, model sensitivity, phase-resolving model, infragravity waves, wave setup, runup wave transformation, Green-Naghdi, model sensitivity, phase-resolving model, infragravity waves, wave setup, runup

Wave flume data from published benchmark experiments were used to extensively evaluate numerical model skill and sensitivity for applying a shock-capturing Green-Naghdi (GN) model to simulate nonlinear wave transformation processes on complex coral reefs. Boussinesq-type models that utilise nonlinear shallow water equations (NSWEs) to represent wave breaking and dissipation hold significant potential for understanding coastal hazards associated with global environmental change and sealevel rise. These fully nonlinear phase-resolving models typically require a threshold condition to switch from dispersive equations to shock-capturing NSWEs in areas of active wave breaking. However, limited information exits regarding how this splitting approach influences the behaviour of different surf-zone processes that contribute to wave runup and inundation on coral reefs. This paper presents a comprehensive analysis of model sensitivity to explore how input parameters that control wave breaking and dissipation influence the behaviour of sea-swell (SS) waves, infragravity (IG) waves, wave setup, runup and solitary waves on coral reefs. Results show that each wave process exhibits unique sensitivity to the free-surface slope threshold (B) that is used to represent areas of active wave breaking by locally switching from the weakly-dispersive GN equations to the shock-capturing NSWEs. However, accurate representation of all wave processes can be achieved if the wave-face steepens to at least 35 degrees (B ≥ 0.7) before breaking is initiated. Results from this research support and encourage the use of nonlinear phase-resolving wave models as tools for academic research, coastal management, coastal engineering and hazard forecasting on atoll and fringing reef environments.

Dear Editor,

We are pleased to submit our article titled Model Skill And Sensitivity For Simulating Wave Processes On Coral Reefs Using A Shock-Capturing Green-Naghdi Solver to be considered for publication in the Journal of Coastal Research. This research paper presents a comprehensive review of model behaviour, skill and sensitivity for simulating nonlinear wave processes on coral reefs in high resolution. Our hope is that this paper opens the door for scientists, managers and engineers to use nonlinear phase-resolving models to examine the impacts of future environmental change to coral reef landforms.

We have prepared all figures in greyscale to suit publication in print and online and have no conflicts of interest.

INTRODUCTION

Coral reefs are characterised by a steep sloping offshore fore-reef that quickly transitions into a shallow and near-horizontal reef flat. This abrupt transition results in energetic wave breaking at the reef edge, where approximately 86% of incident wave energy is dissipated [START_REF] Ferrario | The Effectiveness of Coral Reefs for Coastal Hazard Risk Reduction and Adaptation[END_REF]. Wave energy is further dissipated through friction on the reef flat, with up to 97% reduction of incident wave energy by the shoreline [START_REF] Ferrario | The Effectiveness of Coral Reefs for Coastal Hazard Risk Reduction and Adaptation[END_REF]. The rate of dissipation through breaking and friction are both influenced by reef depth, leading to consistent measurements of tidally modulated wave heights at the shoreline [START_REF] Kench | Wave Processes on Coral Reef Flats: Implications for Reef Geomorphology Using Australian Case Studies[END_REF]. The dependence of shoreline wave height on reef depth means that sea-level rise (SLR) is expected to increase wave energy impacting reef fringed landforms, resulting in widespread flooding and coastal erosion within the 21 st century [START_REF] Storlazzi | Many Atolls May Be Uninhabitable within Decades Due to Climate Change[END_REF].

However, incident wave attenuation on coral reefs is associated with a transfer of energy to free-propagating infragravity waves [START_REF] Péquignet | Energy Transfer between Wind Waves and Low-Frequency Oscillations on a Fringing Reef, Ipan, Guam[END_REF][START_REF] Pomeroy | The Dynamics of Infragravity Wave Transformation over a Fringing Reef[END_REF] and an increase in reef flat water level through wave setup [START_REF] Becker | Water Level Effects on Breaking Wave Setup for Pacific Island Fringing Reefs[END_REF][START_REF] Vetter | Wave Setup over a Pacific Island Fringing Reef[END_REF]. The behaviour of these secondary wave motions on coral reefs is not comprehensively understood and this knowledge gap presents a significant limitation for process based management and engineering on reef coastlines facing SLR. The significance of these secondary wave motions is evident in the fact that all documented cases of wave overtopping on reef landforms identify infragravity waves as the primary process contribution to flooding [START_REF] Cheriton | Observations of Wave Transformation over a Fringing Coral Reef and the Importance of Low-Frequency Waves and Offshore Water Levels to Runup, Overwash, and Coastal Flooding[END_REF][START_REF] Ford | Reef Flat Wave Processes and Excavation Pits: Observations and Implications for Majuro Atoll, Marshall Islands[END_REF][START_REF] Roeber | Destructive Tsunami-Like Wave Generated by Surf Beat over a Coral Reef During Typhoon Haiyan[END_REF][START_REF] Shimozono | Combined Infragravity Wave and Sea-Swell Runup over Fringing Reefs by Super Typhoon Haiyan[END_REF].

Accurate modeling of wave transformation processes, swash dynamics, runup and overtopping on coral reefs landforms is critical for understanding the physical vulnerability of reef coastlines facing SLR and environmental change [START_REF] Reyns | A report to the Asian Development Bank within the framework of the ADB -UNESCO -IHE Knowledge Partnership[END_REF]. However, for models to be an effective hazard management tool, they must accurately represent all the nonlinear wave transformation processes that influence runup and inundation [START_REF] Roeber | Destructive Tsunami-Like Wave Generated by Surf Beat over a Coral Reef During Typhoon Haiyan[END_REF]. Previous attempts to model the impact of SLR on reef landforms (e.g. [START_REF] Storlazzi | Many Atolls May Be Uninhabitable within Decades Due to Climate Change[END_REF][START_REF] Storlazzi | Wave-and Tidally-Driven Flow and Sediment Flux across a Fringing Coral Reef: Southern Molokai, Hawaii[END_REF] rely on phase-averaging techniques that oversimplify energy transfer from incident waves to secondary surf-zone processes and do not simulate high frequency water level motions associated with SS waves, runup or overtopping. Recently, [START_REF] Roeber | Destructive Tsunami-Like Wave Generated by Surf Beat over a Coral Reef During Typhoon Haiyan[END_REF] show that phase-averaging models typically used in coastal hazard assessments do not resolve nonlinear energy transfers in the surf-zone or wave generated flooding. A 'free-surface' or phase-resolving model is required to directly simulate all surfzone processes, runup and overtopping motions.

Coral reef morphology creates an inherent challenge for applying phase-resolving wave models because there is a rapid transition from dispersion dominant processes on the fore-reef slope to flux dominant processes at the reef edge [START_REF] Fang | Fully Nonlinear Modeling Wave Transformation over Fringing Reefs Using Shock-Capturing Boussinesq Model[END_REF]Roeber and Cheung, 2012b). The first challenge is to propagate waves through the shoaling zone in a way that captures the asymmetrical transformation in wave shape and height prior to breaking [START_REF] Huntley | Waves[END_REF]. The second challenge is to resolve wave breaking processes in a way that represents wave shape, dissipation and the generation of secondary wave motions (Roeber and Cheung, 2012b). To overcome these challenges, recent advances have focused on developing fully nonlinear Boussinesq-type solvers that apply 'shock-capturing' techniques to locally represent wave breaking using nonlinear shallow water (NSW) equations (Bonneton et al., 2011a;Bonneton et al., 2011b;[START_REF] Lannes | A New Class of Fully Nonlinear and Weakly Dispersive Green-Naghdi Models for Efficient 2D Simulations[END_REF][START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF]Roeber and Cheung, 2012b;Tissier et al., 2012a). A number of recent studies have shown that shock-capturing Boussinesq-type models can represent incident wave shoaling, breaking and attenuation on coral reefs, while also replicating the generation of infragravity waves and wave setup [START_REF] Beetham | Wave Transformation and Shoreline Water Level on Funafuti Atoll, Tuvalu[END_REF][START_REF] Fang | Fully Nonlinear Modeling Wave Transformation over Fringing Reefs Using Shock-Capturing Boussinesq Model[END_REF]Roeber and Cheung, 2012b;[START_REF] Shimozono | Combined Infragravity Wave and Sea-Swell Runup over Fringing Reefs by Super Typhoon Haiyan[END_REF][START_REF] Su | Boussinesq Modeling of Spatial Variability of Infragravity Waves on Fringing Reefs[END_REF]. These studies are all based on different solutions for the Boussinesq and NSW equations and collectively highlight the capability of phase-resolving models to simulate nonlinear wave processes on complex coral reef environments. Each of these models also applies a threshold condition for wave instability that is used to locally and temporarily switch from dispersive equations to shock-capturing equations in areas of active wave breaking. Complex techniques for switching between source terms include: a momentum gradient threshold (Roeber and Cheung, 2012b); a Froude number threshold [START_REF] Su | Boussinesq Modeling of Spatial Variability of Infragravity Waves on Fringing Reefs[END_REF]; and a wave height to water depth threshold [START_REF] Shimozono | Combined Infragravity Wave and Sea-Swell Runup over Fringing Reefs by Super Typhoon Haiyan[END_REF]. However, the simplest approach that represents how wave breaking is physical triggered is a threshold freesurface slope [START_REF] Beetham | Wave Transformation and Shoreline Water Level on Funafuti Atoll, Tuvalu[END_REF]Bonneton et al., 2011b;[START_REF] Fang | Fully Nonlinear Modeling Wave Transformation over Fringing Reefs Using Shock-Capturing Boussinesq Model[END_REF][START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF]. Regardless of the method, the threshold value chosen to initiate wave breaking has a significant influence on breakpoint location and therefore directly influences the generation and behaviour of secondary wave processes. However, the exact sensitivity of different surfzone processes to wave breaking parameterisation remains unknown.

The aim of this research is to quantify the sensitivity of different surf-zone processes that contribute to wave runup and overtopping on coral reefs to the parameters that influence wave breaking and dissipation in the open source Green-Naghdi model from [START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF].

The Green-Naghdi (GN) equations are a large amplitude solution of the Boussinesq equations which makes them particular suitable to coastal reef environments [START_REF] Fang | Fully Nonlinear Modeling Wave Transformation over Fringing Reefs Using Shock-Capturing Boussinesq Model[END_REF].

With the growing application of shock-capturing Boussinesq-type models on coral reef environments, it is important for users to understand how threshold values chosen to initiate wave breaking influence the full spectrum of wave transformation processes. This paper presents results form a comprehensive sensitivity analysis of how the free-surface slope threshold value (B) influences SS wave dissipation, IG wave behaviour, setup magnitude, runup elevation and solitary wave dynamics on morphologically complex coral reefs. The combined influence from friction is accounted for by testing the sensitivity of B when combined with different friction coefficient (Cf) values. A series of published benchmark scenarios were used to undertake this analysis, including: include irregular wave transformation on a 2D reef [START_REF] Demirbilek | Laboratory Study of Wind Effect on Runup over Fringing Reefs[END_REF]; solitary wave interaction with a 2D reef (Roeber and Cheung, 2012b); and solitary wave interaction with a 3D reef [START_REF] Lynett | Experimental Study of Solitary Wave Evolution over a 3D Shallow Shelf[END_REF].

METHODS

This research uses the shock-capturing Green-Naghdi model developed by [START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF] to simulate laboratory scale scenarios of wave interaction with coral reefs. The Green-Naghdi solution is an attractive option for coral reefs because the weakly-dispersive equations can be applied to relatively large amplitude waves  

Bh

, where a is the characteristic wave amplitude, ho is water depth and Bz is the topographic gradient [START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF]. These characteristics provide a fully nonlinear solution of wave behaviour in shallow water, but are limited to weakly-dispersive scenarios relatively close to the shoreline   , where Lo is the offshore wave length (Bonneton et al., 2011b[START_REF] Lannes | A New Class of Fully Nonlinear and Weakly Dispersive Green-Naghdi Models for Efficient 2D Simulations[END_REF][START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF]. For a description of the governing equations and numerical scheme of the model used in this paper, the reader is referred to [START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF] or the Basilisk website [START_REF] Popinet | A Solver for the Green-Naghdi Equations[END_REF], where the documented source code and a series of examples are freely available.

Wave Breaking

The weakly-dispersive finite-difference solver for wave propagation and shoaling is locally removed in areas of wave breaking, allowing wave shape and dissipation to be solved with a finite-volume solution of the NSW equations [START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF]. Similar approaches to represent wave breaking are becoming increasingly popular in a new generation of shock-capturing GN models (Bonneton et al., 2011a;Bonneton et al., 2011b;[START_REF] Lannes | A New Class of Fully Nonlinear and Weakly Dispersive Green-Naghdi Models for Efficient 2D Simulations[END_REF][START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF]Tissier et al., 2012a;Tissier et al., 2012b). In [START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF], the weakly dispersive GN equations are used unless the local free-surface slope     exceeds a user defined threshold (B), which by default is 1 (45°). At each time-step, on all cells where   > B is true (as calculated by the GN equations), the dispersive function is locally deactivated and the finite volume NSW solution is applied to resolve wave shape and dissipation throughout the breaking process. Therefore, the dispersion equations may be re-applied shoreward of the breakpoint if reef flat water depth can accommodate broken wave heights to a point that the free-surface slope decreases below the B threshold. This is important for accurately propagating reformed wave energy across the reef flat, at which point friction becomes the main form of dissipation. For the simulations presented here, an implicit quadratic bottom friction term [START_REF] Merrifield | Observations and Estimates of Wave-Driven Water Level Extremes at the Marshall Islands[END_REF] was included as:

S u u, f f C  (1)
where Cf is a non-dimensional coefficient that influences the rate of dissipation across the reef flat.

Application to Coral Reefs

Shock-capturing GN models have been used to simulate wave transformation, infragravity waves and setup on laboratory scale [START_REF] Fang | Fully Nonlinear Modeling Wave Transformation over Fringing Reefs Using Shock-Capturing Boussinesq Model[END_REF] and field scale [START_REF] Beetham | Wave Transformation and Shoreline Water Level on Funafuti Atoll, Tuvalu[END_REF] coral reefs. Simulations presented in [START_REF] Beetham | Wave Transformation and Shoreline Water Level on Funafuti Atoll, Tuvalu[END_REF] show that the GN model from [START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF] can accurately represent field scale wave transformation processes. However, the [START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF] model has yet to be evaluated against benchmark data for wave runup and overtopping on coral reefs. Therefore, the [START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF] GN model was used in this research to provide a new and comprehensive evaluation of the model skill and sensitivity for simulation nonlinear wave transformation processes, runup and overtopping processes on a range of coral reef environments. Such benchmarking is required for confidently using the model as a research and coastal hazard management tool.

A visual description of how the shock-capturing GN model represents wave processes on coral reefs is presented in Figure 1, in the context of irregular wave transformation on an idealised reef morphology. This example highlights how the shock-capturing scheme is implemented near the reef edge to produce an asymmetrical wave bore profile on the outer reef flat, characterised by a vertical face and gradual tail. Figure 1 also highlights how incident wave attenuation is associated with setup of mean water level on the reef flat (Figure 1a) and a distinct transfer of energy into the IG wave band (Figure 1d-e). This example shows that maximum wave runup elevation is a function of remnant SS wave energy, IG wave oscillations and wave setup.

Outline of Test Cases

Benchmark test one (BM1) simulated wave transformation processes across a 'typical' Guam fringing reef (Figure 2a) from the University of Michigan (UM) wind wave flume experiment presented by [START_REF] Demirbilek | Laboratory Study of Wind Effect on Runup over Fringing Reefs[END_REF]. The UM experiment was initially undertaken to evaluate the Boussinesq model by Demirbilek and Nwogu (2007), and has subsequently become a standard benchmark for applying wave models to coral reef environments (e.g. [START_REF] Buckley | Evaluation of Nearshore Wave Models in Steep Reef Environments[END_REF]Demirbilek and Nwogu, 2007;[START_REF] Demirbilek | Laboratory Study of Wind Effect on Runup over Fringing Reefs[END_REF][START_REF] Filipot | Spectral Wave Modeling in Fringing Reef Environments[END_REF][START_REF] Nwogu | Infragravity Wave Motions and Runup over Shallow Fringing Reefs[END_REF][START_REF] Sheremet | Modeling of Nonlinear Wave Propagation over Fringing Reefs[END_REF][START_REF] Shimozono | Combined Infragravity Wave and Sea-Swell Runup over Fringing Reefs by Super Typhoon Haiyan[END_REF][START_REF] Su | Parametric Wave-Breaking on Steep Reefs[END_REF][START_REF] Zijlema | Modelling Wave Transformation across a Fringing Reef Using SWASH[END_REF]. Here, the [START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF] GN model was evaluated against fifteen incident scenarios from the UM dataset. The model was used in 1D, with measured water level from wave gauge 1 (WG1) imported as the wave boundary condition. Model outputs for SS wave height (Hss), IG wave height (Hig), wave setup ( ) and maximum runup (Rmax) were compared with flume data to assess model accuracy and sensitivity. A comprehensive sensitivity analysis was undertaken to identify how systematic variations in B and Cf influence each output parameter (Hss, Hig,  and Rmax). Results presented below identify the B and Cf combinations that best represents the range of processes across the surf and swash zone, out of 14 Cf values and 8 B values (112 combinations) for the 15 test conditions (1,680 simulations) Roeber and Cheung (2012) were used to evaluate model results.

The HI reef test simulated a 0.75 m solitary wave breaking on an emerged reef crest, with water plunging into and propagating across a shallow lagoon (Figure 2b). BM2 tested numerical representations of dispersion, nonlinearity, wave breaking, overtopping and wet/dry dynamics (Roeber and Cheung, 2012b). Bathymetry and water level data from the wave flume experiment were downloaded from http://hydraulic.lab.irides.tohoku.ac.jp/app-def/S-102/2014/?page_id=56.

Benchmark test three (BM3) was taken from the National Tsunami Hazard Mitigation

Program benchmarking repository (NTHMP, 2012). BM3 simulated a solitary wave interacting with a complex 3D shelf that has a triangular extension and an emerged conical island (Figure 2c), from the flume experiments by [START_REF] Lynett | Experimental Study of Solitary Wave Evolution over a 3D Shallow Shelf[END_REF]. BM3 used the GN model in 2D (depth averaged) and required solutions for dispersion, shoaling, refraction, diffraction, breaking, overtopping and wet/dry dynamics. The 3D shelf experiment has become a popular Boussinesq model benchmark scenario [START_REF] Fang | An Efficient Shock Capturing Algorithm to the Extended Boussinesq Wave Equations[END_REF][START_REF] Shi | A High-Order Adaptive Time-Stepping Tvd Solver for Boussinesq Modeling of Breaking Waves and Coastal Inundation[END_REF]Roeber and Cheung, 2012a;[START_REF] Yamazaki | NEOWAVE[END_REF] with bathymetry and wave gauge data available online (http://coastal.usc.edu/currents_workshop/problems.html). The same nine combinations of Cf and B from BM2 were tested again in BM3.

Output Data Analysis

For BM1, time-series data for free-surface water level from the UM wave flume and numerically modeled water level outputs at each WG location were analysed using the same (2)

A skill value of one indicates a perfect match between predicted and observed data and a skill value of zero indicates that there is no agreement between the predicted and the observed [START_REF] Lowe | Wave-Driven Circulation of a Coastal Reef-Lagoon System[END_REF]. Model skill was calculated individually for each wave statistic (Hss, Hig, setup, runup) in BM1. Mean skill averaged across these output parameters was then used to understand overall model accuracy and sensitivity for representing wave processes at the shoreline. For BM2 and BM3, model skill was calculated using time-series data for the freesurface water elevation, from numerical and wave flume outputs at each wave gauge location.

RESULTS

This section describes the results from each of the benchmark tests and presents a comprehensive analysis of model behaviour on coral reefs. Collectively, these results identify values for B and Cf that are suitable for accurately simulating nonlinear wave transformation processes, runup and overtopping on coral reef environments.

Benchmark Test 1: Wave Transformation and Runup

Model outputs from BM1 show that the shock-capturing GN equations are capable of replicating flume measurements of SS wave height, IG wave height, setup and runup with maximum model skill between 0.90 and 0.99 (Table 1). This section outlines the sensitivity of each individual wave process to variations in B and Cf and identifies values that provide the most accurate representation of all wave processes at the shoreline, including runup elevation.

Model Sensitivity to the Breaking Slope Threshold (B)

The breaking slope threshold influences the water depth where wave breaking is triggered and determines whether waves break on the reef slope, reef edge or reef flat. The B value therefore influences the rate and mode of incident wave dissipation which has implications for IG wave behaviour and the magnitude of wave setup on the reef flat. Maximum model skill for Hss at the shoreline was 0.90, achieved using B = 0.4 (Figure 3a; Table 1).

However, skill = 0.89 was achieved using B values between 0.4 and 1.2, indicating that a range of values can be used to predict shoreline wave height (Figure 3a). Model outputs for Hss at the shoreline (WG9) were consistently higher using B = 0.4, with a slight decrease in skill and wave height observed when B increased from 0.4 to 0.7 (Figure 3a; Figure 4a-i). However, Hss values were highly variable using these lower B values, with 11% difference in outputs between B = 0.4 and B = 0.7. In comparison, consistent outputs for Hss were achieved using higher B values (0.7-1.2), characterised by 0.9% variation (Table 1). Overall, maximum skill was associated with B = 0.4 but high skill and minimum sensitivity were achieved for Hss using B ≥ 0.7, especially when combined with Cf ≤ 0.008.

Infragravity waves were highly sensitive to B, with a 16.7% difference between minimum and maximum values for Hig at the shoreline (Figure 3b; Table 1). Interestingly, this range decreased dramatically to only 1% variability between Hig outputs using B ≥ 0.7. IG wave height at the shoreline was highest using B = 0.4, with a linear decrease in Hig observed as B increased to 0.7 and consistent outputs between B = 0.7-1.2 (Figure 4a-i). However, maximum model skill (0.98) for shoreline Hig was achieved across a range of B values between 0.4 and 0.8 (Figure 3b). Lower B slopes achieved maximum skill when combined with high friction, compared to higher B values generating maximum skill across a range of low and moderate friction values (Figure 3b). Overall, it was possible to achieve high model skill using all B values (depending on friction) but IG behaviour was highly variable using B ≤ 0.6, compared to consistent Hig predictions with minimal sensitivity associated with B ≥ 0.7 (Figure 3b).

Model outputs for setup magnitude at WG9 were highest using B = 0.4, with a 15.1% decrease in setup observed between B = 0.4 and B = 0.7. Higher B values (≥0.7) produced consistent setup magnitudes, with minimal variation between B = 0.7 and B = 1.2. Low B values over-predicted setup and were associated with a 7.5% decrease in model skill. The increase in setup using low B values explains why Hig and Hss were highest and often over-predicted using B ≤ 0.6. Maximum skill for wave setup (0.99) was achieved using B = 1 but skill above 0.98 was achieved across a wide range of friction values when combined with B ≥ 0.7 (Figure 3b).

Maximum runup elevation is influenced by the magnitude of Hss, Hig and setup at the shoreline. Therefore, model outputs for Rmax reflect the combined sensitivity of all surf-zone processes. Hss, Hig and setup were all largest using B = 0.4, which resulted in significantly overpredicted wave runup outputs (Figure 4j-r) and low skill (Figure 3e), unless unrealistically high friction values were used. Maximum Rmax skill (0.94) was achieved using B = 0.7 (combined with Cf = 0.01), with peak skill for most friction coefficients associated with B = 0.7. Runup skill above 0.92 was achieved across a range of B values between 0.6 and 1.2, with higher friction required to produce accurate runup results when combined with low B values. Runup elevation was 20% lower (on average) and much closer to flume measurements when using B ≥ 0.7.

Model Sensitivity to Friction Coefficient

Quadratic bottom friction has a direct influence on wave attenuation across the model domain and has a significant impact on SS and IG wave transmission across the reef flat.

Maximum model skill for Hss (0.90) was associated with Cf = 0.004, with wave heights typically over-predicted using lower values and under-predicted using higher values (Figure 4a-i). SS waves were reasonably sensitive to Cf, with an average 14.6% decrease in shoreline Hss observed as Cf increased from 0.001 to 0.01 (Table 1). Friction values above 0.01 resulted significantly under-predicted wave heights (Figure 4) and a substantial decrease in skill (Figure 3a). Maximum model skill for Hig was achieved with Cf = 0.006 (0.98), but skill > 0.95 was achieved across a range of Cf values between 0.002 and 0.008 (Figure 3b). IG waves were highly sensitive to Cf, with an average 23.2% decrease in shoreline Hig as Cf increased from 0.001 to 0.01. IG wave transmission across the reef flat was highly influenced by Cf, with low friction resulting in shoreward amplification and high friction resulting in shoreward dissipation. Friction had no influence on wave setup, with only 0.7% variability between outputs and skill > 0.98 achieved by all Cf values (Figure 3c). Wave runup was highly sensitive to Cf, with a 31% decrease in Rmax observed when Cf increased from 0.001 to 0.01 (Table 1).

Rmax was typically over-predicted using Cf ≤ 0.06 and under-predicted using Cf ≥ 0.02, with close predictions of flume runup and high skill (>0.92) achieved using Cf = 0.005 -0.01 (Figure 4j-r). Runup results reflect how SS waves and IG waves respond to friction, with low friction associated with slightly larger SS waves and significantly larger IG wave heights at the shoreline. Runup sensitivity to Cf is primarily associated with how friction influences IG wave transmission across the reef flat.

Maximum Shoreline Skill

Model skill values calculated for each individual wave processes were most consistent when using higher B values (0.7 -1.2) with moderate friction (Cf = 0.002 -0.008). To quantify the accuracy and sensitivity of surf-zone processes, mean skill at WG9 was calculated by averaging skill values for Hss, Hig and setup. The best representation of all wave processes at WG9 (mean skill = 0.95) was achieved using B = 0.8 and Cf = 0.003, with skill > 0.94 associated with all B values between 0.6 and 1.2 combined with Cf between 0.001 and 0.006 (Figure 3d). Wave processes at WG9 were highly sensitive to low B values, with a 14% decrease in the sum amplitude of Hss, Hig and setup between B = 0.4 and B = 0.7. Higher B values produced consistent outputs with less than 1% variation between B = 0.7 and B = 1.2. Output sensitivity to low B values was most pronounced for IG waves (17%) and setup (16%), with SS waves observing slightly less variation (12%). Overall, model skill at WG9 was higher and outputs were more consistent when using B values between 0.7 and 1. Increasing friction from 0.001 to 0.01 had no impact on setup but did force a linear decrease in SS and IG wave height at WG9 (Figure 4). Friction had a major influence on IG wave behaviour on the reef flat, with a 23% decrease in Hig recorded between Cf =0.001 and Cf = 0.01, compared to a 15% decrease in Hss (Table 1). Overall, the most accurate representation of all wave processes at WG9 (skill > 0.94) was achieved using Cf values between 0.001 and 0.006, combined with B ≥ 0.7.

Mean skill at WG9 was combined with runup skill to identify input values for Cf and B that provide the best representation of all shoreline processes that are of key interest to coastal management applications. Maximum shoreline skill was associated with B = 0.7 and Cf = 0.005, with consistent outputs and high skill (>0.93) achieved using B = 0.7-1, combined with Cf = 0.004 -0.007 (Figure 3f). Including Rmax in the assessment of shoreline skill shifted the location of maximum skill to a lower B value and a slightly higher friction value, compared to maximum skill at WG9. However, the location of maximum skill at Cf = 0.005 and B = 0.7 is within the maximum 5% skill contour for each individual wave process (Figure 3f).

Spectral Wave Transformation

Systematic analysis of model skill revealed that the most accurate representation of all wave processes at the shoreline was achieved using B = 0.7 and Cf = 0.005. However, calibrating the model to achieve the best representation of shoreline wave processes does not necessarily confirm that wave processes across the reef platform were accurately represented.

Therefore, wave spectra calculated from numerical model outputs were compared with flume measurements at different wave gauge locations across the reef to further evaluate the accuracy of using B = 0.7 and Cf = 0.005 (Figure 5). Spectra results show a generally strong agreement between the numerical model and the wave flume at the reef slope (WG5) and across the reef flat (WG7-9). These spectral density plots highlight that the GN model adequately represents the dissipation of SS wave energy across the reef and captures the transfer of energy into longperiod IG motions on the reef flat. Further, spectra calculated from GN model outputs show a close match with wave flume spectra across a wide range of incident wave conditions and reef depths that represent high tide, mean tide, low tide and dry reef scenarios (Figure 5).

Benchmark Test 2: Solitary Wave Breaking on a 2D Reef

The HI reef benchmark experiment was simulated using three Cf values (0.001, 0.005, 0.01) combined with three B values (0.4, 0.7, 1), a total of nine model runs. These values encompass the range tested in BM1 and include the combination of B = 0.7 and Cf = 0.005 that achieved maximum skill at the shoreline. This test requires a solitary wave to pass over a dry reef crest (WG10) and is therefore a fundamental benchmark scenario for wave overtopping processes on coral reefs. Skill values reported in this section are based on comparison of freesurface water level between the numerical model and wave flume.

Model Skill and Sensitivity

All combinations of B and Cf were capable of accurately representing wave shape and timing seaward of the reef crest (WG1-8), with mean skill across the first 8 gauges between 0.96 and 0.99 (Figure 6). Variation in skill throughout the shoaling process was attributed to friction, with maximum skill associated with Cf = 0.005 and slightly less skill with higher or lower friction (Figure 6). The value used for B had less influence on offshore skill with B = 0.4 and 0.7 producing near identical wave profiles and skill outputs (Figure 6). B = 1 was associated with slightly less skill between WG1 and WG8 (Figure 6). The wave face became vertical at the upper reef slope (WG9) before breaking on the reef crest (WG10) between t = 12-15 s (Figure 7). Maximum model skill at WG9 was associated with Cf = 0.005, with skill = 0.99 achieved by all B values. Skill at WG9 was slightly lower using Cf = 0.01 and much lower using Cf = 0.001. For all input combinations, model skill decreased between WG9 and WG10 and continued to decrease until WG12, before increasing again at WG13-14 (Figure 6). The highest skill between WG10 and WG12 was consistently achieved using Cf = 0.005, with skill > 0.9 achieved across all wave gauges, regardless of B. Maximum skill in the lagoon was located at WG13 and WG14, with skill > 0.95 achieved using Cf = 0.001 and Cf = 0.005 (Figure 6).

The B value associated with maximum skill at each wave gauge in the lagoon was variable, but mean skill between WG9 and WG14 was slightly higher using B = 0.7 (0.948) compared to B = 0.4 (0.945) or B = 1 (0.947). Overall, mean skill across all wave gauges was maximum using Cf = 0.005 and B = 0.7, which is consistent with outputs from BM1.

Model outputs from BM2 show that the shock-capturing GN equations can replicate the highly nonlinear scenario of at 0.75 m solitary wave breaking on an emerged reef crest.

Using B = 0.7 and Cf = 0.005, skill was calculated separately for the initial forward propagating wave (t = 0-20 s) and the reflected water level surges (t = 20-60 s). Forward wave propagation and shoaling between WG1 and WG8 was simulated with skill > 0.997, consistent with the numerical wave producing a close visual match to time-series outputs form the wave flume (Figure 7). Numerical representation of wave breaking between WG9 and WG10 was characterised by skill > 0.988, confirming the models ability to represent wave shape, amplitude and timing. A significant decrease in forward propagating skill (0.775) was observed at WG11, where the broken wave initially plunged into the lagoon creating a 0.37 m high spike in water level that lasted for 0.4 s in the wave flume. Numerical outputs at WG11 did create a similar surge event but the timing was slightly delayed, the amplitude was higher (0.45 m) and the spike remained at peak amplitude for longer (0.58 s) compared to flume measurements.

Forward propagating model skill increased again between WG12 and WG14, indicating that surge timing, amplitude and shape was reasonably well represented before the wave reflected off the right boundary at x = 87 m. The reflected bore propagated back across the lagoon between t = 20-30 s, with the reflected signal characterised by skill > 0.9 at all lagoon sensors (WG11 -WG14). Less skill (0.87) was observed at WG10 as the surge passed back over the reef crest, but the timing and amplitude of the modeled reflection remained consistent with flume measurements. Flume data show that a series of short period oscillations were created as the reflected wave propagated into deeper water on the upper reef slope. The reflected signal at WG9 was close to flume measurements (skill = 0.99) showing similar timing and amplitude for the short period waves and longer period surge. Skill decreased slightly as the reflected waves propagated towards WG2 (skill = 0.95) but the water level dynamics associated with the reflected wave propagating seaward were accurately captured in the numerical model (Fig. 7).

Benchmark Test 3: Solitary Wave Breaking on a 3D Reef

The 3D shelf with a conical island experiment from [START_REF] Lynett | Experimental Study of Solitary Wave Evolution over a 3D Shallow Shelf[END_REF] is one of the only published benchmark tests that provide freely available data for evaluating 2D phaseresolving wave models. Benchmark test 3 was simulated using the same nine combinations of B and Cf that were used in BM2. BM3 is a popular test for fully nonlinear Boussinesq models because it requires numerical solutions for wave dispersion, shoaling, breaking, diffraction, refraction, convergence and overtopping on a dry shoreline. Snapshots of the model freesurface show how the solitary wave washes over and refracts around the conical cone before surging up the beach-face and overtopping the berm (Figure 8).

Model Skill and Sensitivity

Maximum skill for each combination of B and Cf was recorded near the input boundary (WG1 and WG4) where skill varied between 0.95 and 0.97. Model skill decreased along the y = 0 m transect (WGs 1, 2, 3, 7), with skill at WG2 (in front of the conical cone) between 0.85 and 0.87 depending on input values (Figure 9). Numerical model outputs at WG2 show the timing an shape of the initial wave passing the sensor were close to measurements form the wave flume and suggest that the decrease in model skill is associated with reflections passing over the sensor between t = 20-40 s (Figure 10b). Model outputs directly leeward of the conical cone (WG3) were sensitive to the timing and shape of the initial surges of water passing over and bending around the cone (Figure 10c). Maximum skill at WG3 was 0.88, associated with Cf = 0.001 and B = 0.4, with minimum skill (0.77) associated with Cf = 0.01 and B = 1. Leeward model performance increased further from the cone, with skill at WG7 between 0.86 (Cf = 0.001; B = 1) and 0.90 (Cf = 0.005; B = 0.4). Sensors across the y = 5 m transect (WGs 4, 5, 6, 8) recorded wave processes lateral to the cone and nearly always produced higher skill when compared against gauges on the y = 0 m transect, especially for forward propagating waves (Figure 9). Mean skill across all sensors was highest using Cf = 0.005 combined with B = 0.4 (0.882) or B = 0.7 (0.880). High skill values using B = 0.7 and Cf = 0.005 are consistent with BM1 and BM2 and confirm the use of these values for laboratory scale simulations of wave dynamics on coral reefs.

Velocity

Velocity data were recorded at WG2 and WG3, capturing flow dynamics on seaward and leeward sides of the conical island (Figure 11). In front of the island, WG2 recorded the velocity signal of the steepening solitary wave, with Ux recording a close match to flume measurements using B = 0.7 and Cf = 0.005 (Figure 11a). Numerical outputs at WG2 also recorded a similar Ux signal as reflected water level motions initially surged back over the sensor at t = 7-18 s. Leeward of the island, flume measurements recorded a sharp increase in Ux as the wave surged over the island at t = 7.5 s (Figure 11b). Model results for Ux at WG3 produced the same plunge, followed by a similar trough, between t = 8.5 and 22 s. Modeled Uy data at WG2 and WG3 produced low velocities and did not identify the minor oscillations observed in the flume until t > 22 s (Figure 11d,e). Data from WG10 measured velocity at the upper reef slope away from the island, with numerical outputs for Ux and the more significant variations in Uy showing a reasonable agreement with flume measurements (Figure 11c,f).

Note, skill was not calculated for velocity outputs due to missing sections of data from each of the wave flume sensors.

DISCUSSION

The three benchmark experiments used to evaluate the [START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF] shock capturing GN solver collectively show that the model can skilfully represent surf-zone processes, runup, water level dynamics and overtopping on topographically complex coral reef environments. Of the three benchmark experiments, BM1 provides a comprehensive assessment for applying the model to simulate nonlinear wave transformation processes that contribute to wave driven runup and inundation. BM2 and BM3 provide more detail on how the model represents nonlinear transformations in wave shape and directly test how the GN equations can be applied to simulate wave overtopping scenarios. Physical phenomenon that contribute to wave runup and inundation on atoll islands include: SS waves, IG waves and setup [START_REF] Merrifield | Observations and Estimates of Wave-Driven Water Level Extremes at the Marshall Islands[END_REF].

Field measurements demonstrate how each of these processes respond differently to incident wave conditions and reef depth [START_REF] Becker | Water Level Effects on Breaking Wave Setup for Pacific Island Fringing Reefs[END_REF][START_REF] Ford | Reef Flat Wave Processes and Excavation Pits: Observations and Implications for Majuro Atoll, Marshall Islands[END_REF][START_REF] Kench | Wave Processes on Coral Reef Flats: Implications for Reef Geomorphology Using Australian Case Studies[END_REF][START_REF] Pomeroy | The Dynamics of Infragravity Wave Transformation over a Fringing Reef[END_REF]. Further, these processes influence the behaviour of each other, complicating a conceptual understanding of how environmental conditions influence runup on reef fringed coastlines. Recent developments in numerical modeling (e.g. [START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF]Roeber and Cheung, 2012b;[START_REF] Su | Boussinesq Modeling of Spatial Variability of Infragravity Waves on Fringing Reefs[END_REF] allow these nonlinear interactions to be simulated under controlled boundary conditions, and provide a unique tool for exploring how physical processes will respond to environmental change. However, appropriate and efficient use of these new modeling techniques relies on a comprehensive understanding of how numerical outputs are influenced by parameters that are typically used for model 'tuning' and calibration [START_REF] Oreskes | Verification, Validation, and Confirmation of Numerical-Models in the Earth-Sciences[END_REF]. Systematic analysis from BM1 revealed how SS attenuation, IG transmission, setup and runup each exhibit individual sensitivity to breaking and friction parameters. Therefore, a combined assessment of all parameters is required to achieve a realistic representation of processes that impact the shoreline. In general, outputs were more sensitive to the fiction coefficient (Cf) than to the slope threshold used for locally removing dispersion (B). However, the model was highly sensitive to B = 0.4 (21.8°) because wave breaking was initiated to far offshore. Wave setup was significantly over-predicted using B = 0.4 which meant that IG waves were not sufficiently dissipated by friction processes on the reef flat. Model outputs observed minimal sensitivity and maximum skill using slope thresholds between B = 0.7 (35°) and B = 1.2 (~50°). These findings are consistent with [START_REF] Buckley | Evaluation of Nearshore Wave Models in Steep Reef Environments[END_REF] where a breaking slope of ~55° was found to give the best prediction of UM measurements when using the nonlinear phaseresolving model SWASH (SWASH default is 25°). Using a high breaking slope can be explained physically in this context because of the steep sloping fore-reef on coral reefs typically produces plunging breakers.

Higher friction values were associated with both SS and IG waves being underpredicted at the shoreline. IG waves were highly sensitive to friction, with Hig over-predicted at WG9 using low friction but dramatically under-predicted using high friction. However, skilful predictions for IG wave height were achieved using a Cf range between 0.002 and 0.008.

Model results highlight the role of friction as a primary control for low frequency wave behaviour on coral reefs. Wave setup was consistently well predicted, with minimal sensitivity to friction, unless B = 0.4 was used. Accounting for Rmax added an additional level of complexity when trying to achieve an accurate prediction of all processes at the shoreline.

Maximum runup was highly sensitive to SS and IG dynamics, leading to over-predicted runup with low friction and under-predicted runup using high friction. However, skilful predictions were achieved using moderate friction (0.005 to 0.01) and any B value above 0.6. The best representation of all processes was found using Cf = 0.005 and B = 0.7. These same input values were capable of skilfully replicating solitary wave transformation on complex 2D and 3D reef structures, with a slight decrease in skill observed when using the model in 2DH (two horizontal dimensions). Collectively, the three benchmark scenarios show that the GN model is capable of simulating the range of wave frequencies and nonlinear dynamics associated with wave transformation on coral reef systems.

CONCLUSIONS

A benefit of fully nonlinear free-surface models is that any feedbacks between different frequency wave processes are directly resolved [START_REF] Shimozono | Combined Infragravity Wave and Sea-Swell Runup over Fringing Reefs by Super Typhoon Haiyan[END_REF]. For example, the attenuation of irregular SS waves is associated with the formation of IG motions and setup that in turn influence the transmission of broken and reformed SS waves across the reef flat. These feedbacks are critical for accurately representing the range of processes that promote elevated sea-level and coastal flooding, as recently highlighted by [START_REF] Roeber | Destructive Tsunami-Like Wave Generated by Surf Beat over a Coral Reef During Typhoon Haiyan[END_REF]. However, the interlinked nature of these models also means that the sensitivity of each process needs to be carefully understood before being applied to real world scenarios. Simulations presented here are the first step towards applying the GN solution to study the complex interplay of processes associated with wave transformation on coral reef coastlines. Results suggest that the combination of weakly dispersive Green-Naghdi equations and nonlinear shallow water equations implemented in [START_REF] Popinet | A Quadtree-Adaptive Multigrid Solver for the Serre-Green-Naghdi Equations[END_REF] and other model packages have significant potential for simulating detailed wave processes on reef environments. Consequently, this type of modeling can become a powerful research tool to explore wave dynamics on reef systems in context of future environmental change. wave conditions under high tide (rows 1-2), mean tide (rows 3-4), low tide (rows 5-6) and dry reef (rows 7-8) conditions. Note that the x and y axes are both log scale. Mean skill value across all WG locations is presented inside each subplot.

Tables

  techniques. Simulations ran for 900 s, with mean depth across the platform stabilising after approximately 100 s of wave activity. Data were extracted at 20 Hz (consistent with flume experiment) and trimmed to only include the final 700 s of wave activity when the wave field was developed across the reef. A spectral filter with a 0.25 Hz (4 s) cut-off was used to separate SS and IG wave components. The 0.25 Hz separation is approximately double the peak incident wave period used in the wave basin experiments. Significant wave heights for SS (Hss) and IG (Hig) components were then calculated using a zero down-crossing routine in Matlab. Wave setup ( ) was calculated as the mean free-surface displacement. Wave spectra were calculated from the raw pressure signal using the first 8192 samples from t = 200 s, with an overlapping hamming window and 16 degrees of freedom. Maximum runup (Rmax) was calculated from numerical model outputs using profile data for maximum water level across the reef flat. Rmax values reported inDemirbilek and Nwogu (2007) were used to evaluate numerical model outputs.Model skill was used to quantify how close numerical outputs were to wave flume measurements in all benchmark tests using equation (2), where Pi is the model value, Oi is the wave flume value and vertical bars indicate absolute values.

Figure 1 .

 1 Figure 1. Example of model capability for simulating irregular wave transformation (Hs = 2.6 m, Ts = 10.5 s) on an idealised coral reef (hr = 1 m) with a 1:6 fore-reef slope and 1:14 beach slope. (a) Across reef variations in free-surface water level (wave profile), mean water level (setup) and maximum water level. (b-e) Water level time-series outputs at different locations across the model domain showing the raw signal and the low pass filtered (0.04 Hz) infragravity wave signal.

Figure 2 .

 2 Figure 2. Model bathymetry for each benchmark test showing the location of wave gauge (WG) sensors used to compare numerically modeled outputs with measured wave basin data. (a) Benchmark 1 is the idealised Guam reef from a series of flume experiments undertaken in the University of Michigan wave laboratory[START_REF] Demirbilek | Laboratory Study of Wind Effect on Runup over Fringing Reefs[END_REF]. (b) Benchmark 2 is the Hawaii Reef (HI reef) bathymetry used in the O.H. Hinsdale Wave Flume(Roeber and Cheung, 2012b). (c) Benchmark 3 uses the model in 2D with a complex reef bathymetry that has a steep slope, a triangular shelf and a conical island[START_REF] Lynett | Experimental Study of Solitary Wave Evolution over a 3D Shallow Shelf[END_REF].

Figure 3 .

 3 Figure 3. Model skill contour lines are presented in grey to show how B and Cf values influence the accuracy of wave processes at the shoreline, individually and collectively, from BM1. Shoreline skill at WG9 is presented for (a) SS wave height, (b) IG wave height and (c) setup at the shoreline. (d) Mean skill at WG9 is the average of individual skill values for Hss, Hig and at the shoreline. (e) Rmax skill. (f) Shoreline skill as the average skill from Hss, Hig,  at WG9 and Rmax. Maximum skill is identified with a black marker and the associated Cf and B value is written inside each plot. Black contour lines indicate Cf and B combinations that produce a skill value within 5% of the maximum skill. Note, the x-axis is a log scale.

Figure 4 .

 4 Figure 4. (a-i) Measured (wave flume) outputs for SS wave height (Hss), IG wave height (Hig) and wave setup are compared to numerical model outputs using a selection of different Cf and B values for the 15 UM test scenarios. (j-r) Numerical model outputs for runup elevation are compared to wave flume measurements for the same selection of B and Cf values. The dashed line represents a perfect fit (y = x) and can be used to visually identify how close model outputs are to flume measurements. Skill values written in plots a-i represent mean skill at WG9, with runup skill presented in plots j-r.

Figure 5 .

 5 Figure 5. Wave spectra calculated from the GN model outputs using Cf = 0.005 and B = 0.7 (black) are compared to wave spectra calculated from wave flume data (grey) at the upper reef

Figure 6 .

 6 Figure 6. Across reef variation in free-surface water level skill at the 14 wave gauge locations from BM2, using different B values (represented using different line symbols) and friction values. (a) Cf = 0.001, (b) Cf = 0.05 and (c) Cf = 0.01. (d) Wave flume bathymetry and wave gauge locations are presented to highlight the spatial influence variation in skill across the reef.

Figure 7 .

 7 Figure 7. (a-l) Free-surface water level across the model domain at different time snapshots are presented to show wave propagation, breaking and plunging processes from the GN model (black line), using Cf = 0.005 and B = 0.7, compared with wave flume measurements at each wave gauge location (grey circles). (m-x) Comparison of time-series water level measurements from the GN model (black) and wave flume (grey) using B = 0.7 and Cf = 0.005.

Figure 8 .

 8 Figure 8. Snapshots of free-surface water level at different times during the 3D shelf simulation for BM3, using B = 0.7 and Cf = 0.005.

Figure 9 .

 9 Figure 9. Comparison of model skill for representing free-surface water level at different wave gauge locations, relative to distance from the incident boundary for the y = 0 m (WGs 1, 2, 3, 7), y = 5 m (WGs 4, 5, 6, 8) at y = 10 m (WG9) transect, for each combination of Cf and B.

Figure 10 .

 10 Figure10. Time-series outputs for free-surface water level at each wave gauge location from the GN model using Cf = 0.005 and B = 0.7 (black) are compared to water level measurements from the wave flume (grey).

Figure 11 .

 11 Figure 11. Modeled velocity outputs from the GN model using Cf = 0.005 and B = 0.7 (black) are compared to wave flume measurements (grey) at WG2 (seaward of conical cone), WG3 (leeward of conical cone) and WG10 (lateral to conical cone).

Figure 1

 1 Figure 1

Figure 2

 2 Figure 2

Figure 3

 3 Figure 3

Figure 4

 4 Figure 4

Figure 5

 5 Figure 5

Figure 6

 6 Figure 6

Figure 7

 7 Figure 7

Figure 8

 8 Figure 8

Figure 9

 9 Figure 9

Figure 10

 10 Figure 10

Figure 11

 11 Figure 11

Table 1 .

 1 Summary of numerical model skill and sensitivity when compared to flume measurements of SS wave height, IG wave height, setup and runup at the shoreline. Values assigned to the WG9 column represent the mean skill across all wave process (Hss, Hig and setup) at the shoreline and shoreline skill is the average skill from all wave processes at WG9 and runup. A statistical summary of model skill is first presented, based on all simulations, to provide an indication of how variable outputs are when inputs are modified. To quantify model sensitivity, the percent difference between minimum and maximum values were calculated for all Cf values, all B values and the higher and lower range of Cf and B values. B and Cf values associated with maximum skill are presented as recommended values for predicting each wave process. The range of B and Cf combinations that can be used while maintaining model skill within 5% of maximum skill are shown as the minimum and maximum B and Cf values that make up this range.

		Hss	Hig	Setup	WG9 (Mean) Runup	Shoreline
	Model skill
	Maximum	0.901	0.977	0.988	0.949	0.944	0.936
	Max 5% range	0.856	0.928	0.938	0.902	0.896	0.889
	Mean	0.858	0.887	0.973	0.906	0.798	0.879
	SD	0.057	0.141	0.025	0.065	0.14	0.07
	Minimum	0.668	0.419	0.902	0.692	0.419	0.664
	Output

sensitivity: range between minimum and maximum values

	B (0.4-1.2)	11.9% 16.7% 15.9%	14.8%	20.0%	16.1%
	B (0.4-0.7)	11.1% 16.7% 15.1%	14.3%	20.0%	15.7%
	B (0.7-1.2)	0.9%	1.7%	1.0%	1.2%	1.2%	1.2%
	Cf (0.001-0.01)	14.6% 23.2%	0.7%	12.8%	30.1%	17.2%
	Cf (0.001-0.005)	7.6% 12.2%	0.3%	6.7%	19.4%	9.9%
	Cf (0.006-0.01)	6.1% 10.0%	0.3%	5.5%	10.8%	6.8%
	Recommended B values
	Maximum skill	0.4	0.5	1	0.8	0.7	0.7
	5% range: minimum	0.4	0.5	0.5	0.5	0.6	0.6
	5% range: maximum	1.2	1.2	1.2	1.2	1.2	1.2
	Recommended Cf values
	Maximum skill	0.004	0.006 0.0001	0.003	0.01	0.005
	5% minimum	0.0001	0.002 0.0001	0.0001	0.005	0.002
	5% maximum	0.008	0.008	0.04	0.01	0.01	0.01

ACKNOWLEDGEMENTS

This research was part of a PhD theses completed by Edward Beetham at the School of Environment, University of Auckland (UoA) and was supported by a UoA Doctoral Scholarship.


Samuel Smith
email: ssmith@utas.edu.au

James Venning

Dean Giosio

Paul A Brandner

Bryce W Pearce

Yin Lu Young

Yin Lu

O Tatin

Cloud cavitation behaviour on a

Keywords:

. INTRODUCTION e FSI characteristics of flow over a li ing body can significantly effect the performance of maritime propulsion and control systems. Recent research into the development of composite propellers [ , ] and active control surfaces [ ] has taken place to exploit the ability to passively tailor geometric aspects of the hydrofoil such as skew and pitch based on the loading distribution [ ]. Not only does this self-adaptive behaviour give the ability to design a more energy efficient propeller, but also delay and mitigate the adverse effects of cavitation. One of these effects is the unsteady loading and induced vibration due to the shedding of cloud cavitation.

e effect of unsteady cloud cavitation on the hydroelastic response of hydrofoils has previously been investigated [ , , , , , ], with recently Pearce et al. [ ] showing that the cavity dynamics can influence the FSI response. In addition to the classical shed vortex induced structural response in single phase flow, there is interaction between the development of cavitation on the structural dynamics [ , ]. ese effects are highlighted in experiments by Akcabay et al. [ ] where increased hydrofoil flexibility was seen to increase the cavity length as well as cause a reduction in the cloud cavitation shedding frequency. Further research also shows that flexibility broadens the induced vibration frequency content potentially leading to severe vibration amplification caused by lock-in [ ]. Increased vibrations also occur when the unsteady cavity closure approaches the hydrofoil trailing edge due to high amplitude load fluctuations caused by periodic shedding of sheet-cloud cavitation. In these investigations into the effect of cloud cavitation on hydroelastic response, there is limited discussion on changes in the cavitation pa ern and shedding mechanisms due to the hydroelastic response.

Sheet and cloud cavitation was first extensively studied by Knapp [ ] observing the detachment/shedding of cloud cavitation from a sheet cavity. Since then, several mechanisms have been identified as the primary instability causing periodic shedding depending on the condition. ese included growth of interfacial instabilities such as Kelvin-Helmholtz waves [ , ], re-entrant jet formation [ , , , , , ] and shock propagation [ , , ]. In a recent study on cloud cavitation about a sphere, all three mechanisms have been observed occuring either under varying flow conditions or as a complex coupled mechanism [ ].

To reduce the complexity of the cavitation dynamics, much research into hydrofoil cloud cavitation has focused on D flows to limit D effects and span-wise variations as shown in figure . is is highlighted in time resolved PIV experiments on a D hydrofoil by Foeth et al. [ ] showing significant cavitation stability sensitivity to D flow effects. Span-wise variations are still observed on D hydrofoils where the span-wise cavity length is seen to be proportional to the stream-wise length [ ]. is relationship can result in span-wise cavity lengths that are compatible with the hydrofoil geometry. In these instances, the shedding cloud cavitation exhibits much stronger periodicity than in other conditions [ ].

Force and tip displacement measurements are presented for a nominally rigid stainless steel and flexible composite hydrofoil experiencing cloud cavitation. Synchronised high speed photography is used to analyse the cavitation behaviour and assess the correlation between the cavity dynamics and forces experienced. e aim of this research is to further the understanding of cloud cavitation about a D hydrofoil and how FSI can influence its behaviour.

. EXPERIMENTAL OVERVIEW . Model Hydrofoil Details

Geometry and mechanical properties of the hydrofoil models has been selected based on modelling the static and dynamic fluid-structure interaction typical of propellers and control surfaces. e chosen geometry was a symmetric (unswept) trapezoidal planform of mm span with a mm tip and mm root chord, providing an aspect ratio = . . e chord length was chosen to be compatible with the mounting to the water tunnel test section and to achieve a chord based Reynolds number Re c = 0.7 × 10 6 . e unswept geometry, in conjunction with a span-wise alignment of the fibre orientation, was intentionally chosen to principally consider bending deformation only of the flexible hydrofoil. A modified NACA section profile with a thicker trailing edge was selected for improved manufacture of the flexible composite model (see Zarruk et al [ ] for further details).

e flexible (composite) model was manufactured as a carbon/glass-epoxy hybrid structure consisting of a poly-olefin scaffold core, T unidirectional carbon fibre and biaxial E-glass fabric used as the key structural components with an outermost fine E-glass basket weave layer to aid surface finish. A full lay-up sequence and construction procedure is detailed in [ ] where the composite hydrofoil model used in the present study is termed the CFRP hydrofoil.

e rigid (stainless steel) model was machined from a Type stainless steel billet with both models manufactured to ± . mm surface tolerance and . µm surface finish.

e response spectrum of both hydrofoils was determined from both impact tests and hydrofoil loading spectra [ ] with results summarized in table . First mode natural frequencies were obtained in air at Hz and Hz, and in water at and Hz for the stainless steel (rigid) and composite (flexible) models respectively.

Hydrofoil Dynamic Properties

Hydrofoil Rigid Flexible First bending mode in air (Hz) First bending mode in water (Hz)

. Experimental Setup

Measurements were carried out in the Cavitation Research Laboratory (CRL) water tunnel at the Australian Maritime College. e tunnel test section is . m square by . m long in which the operating velocity and pressure ranges are to m/s and to kPa absolute respectively. e tunnel volume is 365 m 3 with demineralised water (conductivity of order 1 µS/cm). e test section velocity is measured from one of two (high and low range) Siemens Sitransp differential pressure transducers models MF -DA -AB -Z and MF -FA -AB -Z (measuring the calibrated contraction differential pressure) with estimated precisions of . and .

m/s respectively. A detailed description of the facility is given in [ ]. As shown in figure , two profiled plates are used to clamp the model within a housing that is a ached to a -component force balance. e hydrofoil, located at the mid length of the test section, extends vertically into the flow through a mm diameter penetration in the ceiling. e penetration is made fair (to 50 µm) using a disk mounted, in this case, on the measurement side of the balance. e fairing disk has a . mm radial clearance to avoid interference with the force measurement.

. Experimental Techniques

Data was obtained for a cavitation number of . and at a velocity of . m/s which corresponds to a Reynolds number (based on mean chord length) of 0.7 × 10 6 with the models at an incidence of 6°.

e cavitation number is defined as

σ = 2(p -p v )/ρU 2
∞ and Reynolds number as Re = U ∞ c/ν, where p is the static pressure at the test section centreline, p v is the vapour pressure, ρ is the water density, U ∞ is the test section velocity, c is the mean chord and ν is the kinematic viscosity of the water. Of the total load vector measured, only the time-varying component of the normal force is presented.

e force balance was calibrated by a least squares fit between a basis vector loading cycle and the outputs giving a 6 × 6 matrix. An estimated precision on all components is less than 0.1%.

e cavitation behaviour was recorded using high speed photography with a HighSpeedStar (LaVision, Germany) mounted on the side of the test section.

e camera was outfi ed with a Nikkor f/ .

mm lens and setup with a magnification factor of . px/mm. High speed images were recorded with a spatial resolution of 1024 × 1024 at , Hz for the rigid hydrofoil where the flexible foil was recorded at , Hz due to data acquisition limits. A previous study [ ] has shown the stainless steel model to be nominally rigid with a maximum tip deflection of less than % of the mean chord compared with % for the composite model.

e force data for the stainless steel model was found to be nominally invariant with Reynolds number for α ≤ 6°. On this basis, tip bending displacement, δ tip , was only recorded for the flexible model with measurements achieved by tracking . mm diameter white dots on the tip of the hydrofoil. Further information on the technique used in similar experiments can be found in [ ]. is was achieved using a HighSpeedStar high speed camera mounted on the bo om of the test section. e camera was outfi ed with a Nikkor f/ .

mm lens where images had a magnification factor of . px/mm. Images were recorded at , Hz with a spatial resolution of 512 × 1024.

e high speed photography was synchronized with the force measurement acquisition by simultaneous triggering from a BNC Model Pulse Generator. Force and tunnel flow data were sampled at , Hz and , Hz, respectively.

. RESULTS AND DISCUSSION

. Force and Tip Displacement Time Series e measured normal force, X, is presented as a dimensionless coefficient, C X = 2X/(ρAU 2 ∞ ), where A denotes the planform area. e frequency content of the X force experienced by either hydrofoil is shown in the power spectral density in figure . Both hydrofoils exhibit a common primary frequency at approximately Hz with secondary frequencies at and Hz for the rigid and flexible hydrofoils, respectively . e common fluctuation at Hz is linked to periodic shedding of cloud cavitation from mid-span as made evident in spacetime plots discussed later.

ese plots also reveal another shedding mechanism towards the tip of the rigid hydrofoil tip /c, and pixel intensity, I, for both the rigid and flexible hydrofoils. Pixel intensity was taken along the span, s, at 0.24s, 0.47s and 0.77s from the root and at % of the local chord for both hydrofoils. Frames taken from the synchronised high speed video show one full shedding cycle at the dominant frequency with the first and last frames coinciding with the maximum force of the primary frequency component. For both hydrofoils, the primary frequency corresponds to a full length cavity being shed. e maximum force occurs just a er the a ached cavity breaks off and a new cavity has just formed. e minimum force coincides with maximum cavity length. Analysing pixel intensity of the rigid hydrofoil we see a rise, plateau and fall at the middle pixel (0.47s) during the shedding cycle shown. e top pixel (0.24s) shows similar behaviour and duration but slightly trails the middle pixel intensity in time.

δ ′ tip /c f g h i j (f) (g) (h) (i) (j)
e bo om pixel (0.77s) exhibits peaks for a much shorter duration compared to the other pixels and primarily occurs when the top and middle pixels are low.

ese pixel intensity traits suggests there are two shedding modes at the top-middle and the bo om that alternate over time. e flexible hydrofoil shows similar trends but due to low temporal resolution, definitive conclusions cannot be made at this stage.

Comparing the force signals, the rigid hydrofoil exhibits a slightly lower mean C X to that of the flexible hydrofoil, . and .

, respectively. However, the rigid hydrofoil exhibits more unsteadiness with the flexible hydrofoil appearing to almost dampen some of the lower amplitude fluctuations. is is reflected in the C X RMS values of . and .

for the rigid and flexible hydrofoil, respectively. e unsteady displacement of the flexible hydrofoil is seen to vary significantly over time having a strong correlation with C X as expected.

. Shedding Mechanisms and Cavity Dynamics rough analysis of the high speed videos, it is evident that the primary shedding mechanism is the classical re-entrant jet. A typical shedding process can be seen in the space-time plot of the rigid hydrofoil (figure ) generated from a line of pixels extracted at a position mm along the span (i.e. at 0.33s). Once the cavity forms, it initially grows at a constant velocity during stage . At a certain point (t ≈ 0.01), the re-entrant jet starts to propagate forward as indicated by a second curve forming in the cavity. At the same instant, the cavity shi s into its second growth phase with a reduced cavity growth speed. As the re-entrant jet approaches the cavity detachment, it starts interacting with the upper surface of the cavity (t ≈ 0.2). is is indicated by the white streaks from the secondary curve of the re-entrant jet in figure . Shortly a er, the jet reaches the cavity detachment, breaking off the a ached cavity, forming a cavitation cloud that is then advected downstream. Following cavity break-off, a new cavity forms soon a er and the cycle starts again.

At the flow conditions examined here (Re c = 0.7 × 10 6 , σ = 0.8 and α = 6°), the NACA hydrofoil forms a rel- atively thin cavity, resulting in a thin jet that initially (see stage growth in figure ) has insufficient momentum to break through to the cavity detachment point due to friction of the adjacent layers [ , ]. A sufficiently thin cavity may also have significant interactions between the upper and lower interfaces of the cavity as surface perturbations become predominant leading to small-scale vapour structures being shed instead of a large-scale cloud [ ]. is can be seen in chord-wise space-time plots (figure ) of several shedding cycles where there is significant variation between each cycle. e interaction of surface perturbations manifest as a rough opaque surface, seen at the top of figure , compared the transparent region shown in figure . Comparison of the rigid and flexible hydrofoils sees that shedding variations over time exist for both hydrofoils but aren't as severe for the flexible hydrofoil.

e span-wise space-time plots taken mm upstream of mid-chord for both hydrofoils (figure ) illustrates how the cloud cavitation varies along the span over time. It is observed that there is significant span-wise variation for both hydrofoils with no uniform (i.e. across the whole span) shedding observed.

is complex cavitation behaviour is due to the interaction of multiple effects. is includes spanwise flow disparity over the hydrofoil due to the tapered and swept geometry causing changes in re-entrant jet direction. Additionally, the nature of the vertical mounted hydrofoil results in cavitation number gradient along the span as well as buoyancy effects on the cavities.

A level of consistency is seen in the periodicity of shedding, but the shedding behaviour of each event varies with clear re-entrant jet observed in some, but not in others. Analysis of the rigid hydrofoil reveals the existence of two shedding modes along the span at frequencies of approximately Hz and Hz for the upper and lower parts of the hydrofoil, respectively, calculated from space-time plots. e flexible hydrofoil also shows signs of two shedding modes, Hz and Hz, with two crescent cut-outs along the span-wise length of the cavity observed, similar to the rigid hydrofoil. Non-dimensionalizing these frequencies using a cavity length based Strouhal number, St = f L c /U ∞ , where L c is the maxi-mum cavity length, the rigid frequencies equate to . and . for the upper and lower parts, respectively, compared to . and . for the flexible hydrofoil.

Comparing the span-wise space-time plots, the flexible hydrofoils secondary shedding mode at the bo om doesn't appear to have as strong a periodicity as the rigid hydrofoil. is could be due to the force induced tip displacement imposed by the larger shedding cavity from the top interfering with the lower shedding physics. ese observations coincide with frequencies present in the X force spectrum (figure ) supporting the strong correlation between cavity dynamics and forces.

. CONCLUSIONS

Preliminary results are presented for the effect of FSI on cloud cavitation about a hydrofoil.

e normal force, tip bending displacement and cavitation behaviour were compared for a nominally rigid stainless steel and flexible composite NACA hydrofoil at Re c = 0.7 × 10 6 , σ = 0.8 and α = 6°. e bending deformations of the flexible hydrofoil were seen to dampen some of the higher frequency fluctuations in the normal force measurements while showing a strong correlation between tip displacement and normal force.

A re-entrant jet was identified as the primary shedding mechanism showing changes in growth and jet speed at various stages in the shedding cycle. Due to the thin cavity, surface perturbations were seen to have significant interaction, sometimes resulting in small-scale vapour pockets being shed instead of large-scale cloud cavitation. e cavitation behaviour is observed to be highly complex due to the D nature of the flow leading to significant span-wise flow disparity. Both hydrofoils exhibited fairly consistent periodic shedding but varying behaviour between each event. Two shedding modes appeared to form along the span with either hydrofoil showing two curved regions in the cavity trailing edge typical of a re-entrant jet. is results in two shedding frequencies for either hydrofoil with both having a primary frequency of Hz and secondary frequencies of Hz and Hz for the rigid and flexible hydrofoil, respectively. ese differences are a ributed to force induced tip bending displacements affecting the cavity dynamics due to changes in the flow field. is is supported by the fact that the observed shedding frequencies matched those present in the X force spectrum.

e fluid-structure interaction phenomena observed for a flexible D hydrofoil experiencing cavitation is highly complex. is complicated behaviour makes it difficult to predict the performance of real world applications such as composite propellers. Further investigation into the phenomena in a wider range of conditions will allow more detailed and accurate predictions, permi ing improved designs of control surfaces and marine propulsors.

Figure .

 . Figure . Cavitation about a NACA stainless steel hydrofoil at σ = 0.8, Re c = 0.7 × 10 6 and α = 6°.

Figure .

 . Figure . Hydrofoil model assembly showing an exploded view of the clamping housing arrangement allowing continuity of the reinforcing fibres for the CFRP models.

Figure .

 . Figure . Narrowband X force power spectral density (PSD) for the rigid and flexible hydrofoils at σ = 0.8, Re c = 0.7 × 10 6 and α = 6°.

Figure .

 . Figure .Sample time series of normal force coefficient with selected images from high-speed imaging for the rigid (top) and flexible hydrofoil (bo om). Simultaneous normalized unsteady tip bending displacement (dot-dashed line) is also presented for the flexible hydrofoil. e curves of pixel intensity are taken at % of the chord for span-wise locations of 0.24s, 0.47s and 0.77s from the root. Data was taken at σ = 0.8, Re c = 0.7 × 10 6 and α = 6°.

Figure

  Figure shows a short time series of C X , the unsteady tip bending displacement to chord ratio, δ ′ tip /c, and pixel

Figure .

 . Figure . Space-time plot of a single shedding cycle of the rigid hydrofoil . showing the key components of the shedding cycle at σ = 0.8, Re c = 0.7 × 10 6 and α = 6°. e flow direction from top to bo om.

Figure .

 . Figure . Chord-wise space-time plots from high speed images taken mm along the span for the rigid (top) and flexible hydrofoil (bo om) at σ = 0.8, Re c = 0.7 × 10 6 and α = 6°. e flow direction is top to bo om.

Figure .

 . Figure . Span-wise space-time plots from high speed images mm upstream of the mid-chord for the rigid (top) and flexible hydrofoil (bo om) at σ = 0.8, Re c = 0.7 × 10 6 and α = 6°. e flow direction is le to right.

Table .

 . Hydrofoil natural frequencies obtained from impact tests and loading spectra for in-air and in-water frequencies, respectively[ ].

  [ ] C. E. Brennen. e dynamic balances of dissolved air and heat in natural cavity flows. Journal of Fluid Mechanics, ( ): -, . [ ] F. Avellan, P. Dupont, and I. L. Ryhming. Generation mechanism and dynamics of cavitation vortices downstream of a fixed leading edge cavity. In th Symposium on Naval Hydrodynamics, pages -, . [ ] S.P. Furness and S. P. Hu on. Experimental and theoretical studies of two-dimensional fixed-type cavities. B. Stutz and J. L. Reboud. Experiments on unsteady cavitation. Experiments in Fluids, ( ): -, . [ ] M. Callenaere, J. Franc, J. Michel, and M. Riondet. e cavitation instability induced by the development of a re-entrant jet. Journal of Fluid Mechanics, : -, . [ ] K. R. Laberteaux and S. L. Ceccio. Partial cavity flows. part . cavities forming on models without spanwise variation. Journal of Fluid Mechanics, : -, . [ ] J. K. Jakobsen. On the mechanism of head breakdown in cavitating inducers. Journal of Basic Engineering, ( ): -, . [ ] G. E. Reisman, Y. C. Wang, and C. E. Brennen. Observations of shock waves in cloud cavitation. Journal of Fluid Mechanics, : -, . [ ] H. Ganesh, S. A. Mäkiharju, and S. L. Ceccio. Bubbly shock propagation as a mechanism for sheet-to-cloud transition of partial cavities. Journal of Fluid Mechanics, : -, . [ ] KL De Graaf, PA Brandner, and BW Pearce. Spectral content of cloud cavitation about a sphere. Journal of Fluid Mechanics, , . [ ] E. J. Foeth, C. W. H. Van Doorne, T. Van Terwisga, and B. Wieneke. Time resolved piv and flow visualization of d sheet cavitation. Experiments in Fluids, ( ): -, . [ ] Y. Kawanami, H. Kato, and H. Yamaguchi. reedimensional characteristics of the cavities formed on a two-dimensional hydrofoil. In ird International Symposium on Cavitation, volume , pages -. Laboratoire des Ecoulements Géophysiques et Industriels, Grenoble, France, . [ ] G. A. Zarruk, P. A. Brandner, B. W. Pearce, and A. W. Phillips. Experimental study of the steady fluid-structure interaction of flexible hydrofoils. Journal of Fluids and Structures, : -, . [ ] P. A. Brandner, Y. Lecoffre, and G. J. Walker. Design considerations in the development of a modern cavitation tunnel. In th Australasian Fluid Mechanics Conference, pages -, . [ ] Y. L. Young, N. Garg, P. A. Brandner, B. W. Pearce, D. Butler, D. Clarke, and A. W. Phillips. Load-dependent bendtwist coupling effects on the steady-state hydroelastic response of composite hydrofoils. Composite Structures, . [ ] P. F. Pelz, T. Keil, and T. F. Groß. e transition from sheet to cloud cavitation. Journal of Fluid Mechanics, : -, .

	Journal of Fluids Engineering, page	,	.
	[ ] Q. Le, J. Franc, and J. Michel. Partial cavities: Global
	behaviour and mean pressure distribution. Journal of
	Fluids Engineering,	: -,	.
	[ ] Y. Kawanami, H. Kato, H. Yamaguchi, M. Tanimura, and
	Y. Tagaya. Mechanism and control of cloud cavitation.
	Journal of Fluids Engineering,	( ): -,	.
	[ ]

ACKNOWLEDGEMENTS

is project was supported by the University of Tasmania, the Defence Science and Technology Group, US Office of Naval Research (Dr.Ki-Han Kim, Program Officer) and ONR Global (Dr.Woei-Min Lin) through NICOP S&T Grant no. N ---. e authors would like to acknowledge the assistance of Mr Steven Kent and Mr Robert Wrigley from the Australian Maritime College for their essential help with se ing up and carrying out the experiments.