Distant-Source Swells Cause Coastal Inundation on Fiji’s Coral Coast

Distant-source swells are known to regularly inundate low-lying Pacific Island communities. Here we examine extreme total water level (TWL) and inundation driven by a distant-source swell on Fiji’s Coral Coast using observations and a phase-resolving wave model (XBeach). The objective of this study is to increase understanding of swell-driven hazards in fringing reef environments to identify the contribution of wave setup and infragravity waves to extreme TWL and to investigate coastal flooding during present and future sea levels. The maximum TWL near the shore was caused by compounding mechanisms, where tides, wave setup, infragravity waves, and waves in the sea swell frequencies contributed to the TWL. Waves and wave setup on the reef were modulated by offshore wave heights and tides with increased setup during low tide and increased wave heights during high tide. Numerical simulations were able to reproduce the mechanisms contributing to the extreme TWL and allowed an estimation of the inundation extent. Simulations of the same swell under the RCP8.5 sea-level rise scenario suggest the area of inundation would increase by 97% by 2100. A comparison between the numerical model, a multiple linear regression model, and two commonly used parametric models reveals that both XBeach and the linear regression model are better suited to reproduce the nearshore wave setup and TWL than the empirical equations. The results highlight the need for customized, site-specific coastal hazard assessments and inundation forecast systems in the South Pacific.


INTRODUCTION
Coastal inundation can cause shoreline changes (e.g., erosion), destruction of infrastructure, contamination of freshwater reserves and agricultural land, and even loss of lives (Hoeke et al., 2013;Storlazzi et al., 2018;Smith and Juria, 2019). It is, therefore, critical for decision makers to have access to detailed coastal hazard and risk information to adequately manage coastal areas. In the tropical Pacific, coastal inundation is commonly associated with unusually high astronomical tides and storm surges from tropical cyclones (McInnes et al., 2014). The risk of inundation from distant-source swells is less well studied (Hoeke et al., 2013). Consequently, regional meteorological services often fail to warn communities of imminent swell-driven inundation (Bosserelle et al., 2015a;Winter et al., 2020). Hoeke et al. (2013) investigated a series of inundation events in low-lying island communities in the western Pacific that were caused by swell waves generated by extra-tropical storms in the north Pacific. Although the flooding caused severe impacts to five Pacific Island nations, tide gauges failed to record evidence of extreme sea levels. The authors argue that tide gauges and satellite altimetry often fail to capture wave contributions to total water level (TWL) and, consequently, miss many inundation events. A recent study projects that many atolls will become uninhabitable over the coming decades due to wave-driven inundation exacerbated by increased sea levels (Storlazzi et al., 2018). Inundation risk assessments and future early warning systems should, therefore, consider the wave contributions to extreme TWLs (Hoeke et al., 2013).
There are two major wave-driven processes that contribute to the TWL on reef-lined coasts. First, dissipation of swells due to breaking causes wave setup (η wave , elevation of the mean still water surface) from the excess momentum flux in the water column (Becker et al., 2014). Second, swell dissipation over coral reefs generates infragravity (IG) waves (bore-like uprush of individual waves on the beach) (e.g., Baldock, 2012;Pomeroy et al., 2012) that significantly contribute to the total wave runup (e.g., Van Dongeren et al., 2013). η wave has been found to reach up to one third of the incident wave heights H s,inc along coastlines of tropical and subtropical islands (Munk and Sargent, 1948;Tait, 1972;Vetter et al., 2010). Both IG waves and η wave have been shown to cause morphological changes (e.g., Kench et al., 2017) and coastal inundation (e.g., Ford et al., 2018) along fringing reef coastlines.
Recent work has investigated the importance of wave contributions to flooding at inter-annual to decadal time scales in the context of sea-level rise (Rueda et al., 2017;Vitousek et al., 2017;Storlazzi et al., 2018;Vousdoukas et al., 2018). Vitousek et al. (2017) investigated projected changes in coastal flooding on a global scale using sea-level projections combined with wave, tide, and storm surge models. The authors estimate a doubling of the frequency of extreme water levels in the tropics by no later than 2050. Vousdoukas et al. (2018) investigated global projections of extreme sea levels based on a probabilistic approach and estimate that, for a large part of the tropics, once-in-a-century extreme sea levels would become annual by 2050. To provide a large-scale estimate of potential future extreme TWLs, these global studies are generally based on empirical equations to calculate η wave and neglect wave run-up. For example, Vousdoukas et al. (2018) calculated η wave using a generic approximation, where η wave is 2% of H s,inc , and Vitousek et al. (2017) calculated η wave based on the 2% exceedance run-up on dissipative beaches (Stockdon et al., 2006). Although these global assessments serve as useful first-order inundation hazard estimates, they neglect wave run-up due to incident IG waves and may significantly misrepresent extreme TWL and inundation hazard. Furthermore, these studies drew their wave information from coarse global wave simulations and neglected regional features, such as nearshore bathymetry, geomorphology, or local winds. To accurately quantify the wave contribution to extreme TWL at a particular coastline, hydrodynamic models at high spatial and temporal resolutions or observational methodologies should be applied (Hoeke et al., 2011(Hoeke et al., , 2013McInnes et al., 2014;Winter et al., 2020).
In the present study, we examine extreme TWL driven by a distant-source swell in May 2018 on Fiji's Coral Coast (Figure 1) using observations from an offshore wave buoy and a pressure sensor on the reef flat. The objective of this study is to increase our understanding of swell-driven hazards in fringing reef environments and to identify the contribution of η wave and IG waves to extreme TWLs and coastal inundation on reef-lined islands. The inundation extent under present and future sea levels is simulated using a high-resolution numerical model (XBeach; Roelvink et al., 2009) to demonstrate how increased sea levels can modify the relative contribution of the different components. Finally, we compare the XBeach simulations with two empirical equations and a multiple linear regression model to investigate their suitability to reproduce the extreme TWL.

Study Site
Fiji's Coral Coast in the south of the main island of Viti Levu is a major tourist area critical for the economy and development of the country (Harrison and Prasad, 2013). Although the area is known to be frequently flooded by distant-source swells, these events often remain undocumented and have received limited attention in the peer-reviewed literature (e.g., Hoeke et al., 2013;McInnes et al., 2014;Bosserelle et al., 2015a).
Viti Levu features a semi-diurnal tide with 1.2 m amplitude. The wave climate is dominated by long-period southwesterly swells generated by storms in the Southern Ocean and Tasman Sea and by shorter period east-southeasterly waves generated by the trade winds. The annual mean significant wave height (H s ) is 2.15 m with a 99th percentile H s of 3.93 m; 90% of times large wave heights (i.e., waves exceeding the 99th percentile) are generated by Southern Ocean swells. The annual mean peak period (T P ) is 13 s (Bosserelle et al., 2015b).

Observations
Offshore wave conditions off Cuvu on Fiji's Coral Coast (Figure 1) were recorded using a Sofar Spotter wave buoy deployed at 200 m water depth. The wave buoy measured the full directional wave spectrum at a frequency resolution of 0.0098 Hz and a frequency range of 0.033-1 Hz (i.e., wave periods between 1 and 30 s) at a sampling rate of 2.5 Hz. Half-hourly bulk statistics (i.e., H s , mean period (T m ), T P and directions θ as well as directional spreads) were automatically calculated onboard and transmitted in near real-time via an iridium satellite connection while the full spectral data was stored internally and downloaded after retrieval.
The nearshore TWL was recorded using an RBR Duo pressure sensor (Figure 1) programmed to continuously sample at 1 Hz. Water levels were corrected for ambient atmospheric pressure (Fofonoff and Millard, 1983) using data from a nearby automated weather station (Figure 1 and Supplementary Figure S1). η wave and IG wave heights were calculated from the water-level time series. η wave was calculated as the difference between the predicted tide and the 1-h mean observed water level. The water level time series was decomposed into short and IG waves using a continuous wavelet transform with a frequency range for short waves in the sea and swell frequency bands between f min = 0.033 and f max = 0.6 Hz (wave periods between 1.67 and 30 s) and IG waves within lower frequencies between f min = 0.0036 and f max = 0.033 Hz (wave periods between 30 s and ∼5 min). Significant wave heights of short (H s,SW ) and IG waves H s,IG were calculated from the partially integrated power spectrum: (1) Tidal harmonics were derived from past (December 17, 2013-September 30, 2014) water-level measurements at a nearby tide gauge in Maui Bay (Figure 1) using the Unified Tidal Analysis and Prediction (UTide) software (Codiga, 2011).

Multiple Linear Regression Model
Statistical models are often applied as computationally efficient alternatives to numerical models (e.g., Camus et al., 2014;Cid et al., 2017). Here, we investigated the relationship between observed incident H s and tide level (η tide ) (predictor variables) and H s,SW , H s,IG , η wave , and TWL (response variables) on the reef by performing a multiple linear regression analysis: where Y is the response variable (H s,SW , H s,IG , η wave , or TWL) and X 1 and X 2 are the predictor variables (offshore H s and η tide , respectively). a 1 , β 1 , . . . , β 3 are the regression coefficients.

Numerical Model (XBeach)
We simulated the swell using the XBeach model (Roelvink et al., 2009) in surfbeat mode. XBeach is a coupled wave and circulation model that resolves the complete wave transformation (e.g., swash hydrodynamics) including wave-group forcing and resulting IG waves through the non-linear shallow water equations. The model combines the time-varying wave-action balance with a dissipation model to derive the wave-group forcing (Roelvink et al., 2009). Surface rollers inducing a shoreward shift in wave forcing are resolved using a roller model. Long-and cross-shore currents as well as wave-current interaction in the wave boundary layer and resulting changes in wave-averaged bed shear stress are also accounted for. Despite being originally developed for mild sloping beaches, XBeach has since been successfully adapted and applied to fringing reef environments (Pomeroy et al., 2012;Van Dongeren et al., 2013;Buckley et al., 2014;Quataert et al., 2015;Rueda et al., 2019). This has been achieved by including a bottom friction dissipation term in the short-wave energy balance and by increasing the bottom friction coefficient for flow cf , which is associated with currents and IG wave action, and the incident short wave dissipation fw (Van Dongeren et al., 2013;Quataert et al., 2015). Furthermore, XBeach was recently improved to better predict the groupiness of the short waves and the resulting IG waves by implementing the "single-dir" option and by improving the propagation scheme for the wave-action balance (Roelvink et al., 2018). In the present study, we apply this new feature to our simulations. We simulated the time span of May 25-29, 2018, and stored the free surface TWL (zs), maximum TWL (maxzs), and Hrms wave height.
The wave model was forced at the boundary with hourly parametric wave conditions (H s , T P , θ P , and directional spread) measured by the wave buoy. η tide was calculated from historic tidal records at Maui Bay (see section "Study Site"). Meteorological sea level anomalies (e.g., wind setup and inverted barometric effect) were not taken into account in the numerical model to isolate wave-driven effects on the reef.
Bathymetry and topography data were obtained from a range of different surveys and data sets: Bathymetry data in the deep water (i.e., off the reef edge) were collected using boat-mounted multi-and single-beam sensors. The bathymetry on the reef platform (shallower than ∼10-15 m) and the topography data were derived from satellite imagery. The dissipation of waves across a reef is extremely sensitive to the geomorphology (i.e., reef slope, water depth, reef width, or bed friction) (Rueda et al., 2019). Initial attempts to validate the model failed due to inaccuracies in the satellite-derived bathymetry on the reef flat. A detailed real-time kinematic (RTK) positioning survey on the reef and along the beach allowed us to correct the satellite data. All data were tied to ground control points, which were referenced to a primary benchmark in Maui Bay, which, in turn, was linked to a tide gauge in Maui Bay. All data (including the water levels in this study) are referenced to the local mean sea level.
All data were collated and interpolated onto a 10-m resolution grid (940 by 515 grid points; Figure 1). As η wave , wave breaking, and other effects, such as the generation of IG waves, is highly dependent on the reef slope, it can be assumed that the response of wave-driven inundation to sea-level rise is non-linear. We, therefore, simulated the May 2018 swell with current sea levels derived from the Maui Bay tide gauge and with increased water levels (0.84 m), corresponding to the likely global mean sea-level rise under the representative concentration pathway (RCP) 8.5 (Portner et al., 2019).
A multivariate sensitivity analysis was performed to find ideal values for the slope of the breaking wave front in the roller model (β), the breaker parameter (γ) (Roelvink, 1993), the bottom friction coefficient for flow cf , and the short wave friction coefficient fw . Previous studies (Lowe et al., 2007;Van Dongeren et al., 2013) indicate that the short wave friction coefficient and bottom friction coefficient for swell should be considered independently. Other parameters were set to default (e.g., Roelvink et al., 2009). The sensitivity analysis was performed as follows: (1)  (3) For each simulation, the root mean square error (RMSE) was calculated forη, H s,SW , and H s,IG . The lowest RMSE averaged acrossη, H s,SW , and H s,IG was found for γ = 0.6, β = 0.3, fw = 0.3, and cf = 0.005, and these values were applied to the simulation of the entire swell. Although the fw value agreed with previous findings, the cf value was lower than expected when comparing to similar studies (e.g., Péquignet et al., 2011). It is well documented that small changes in reef composition can have large effects on the bottom roughness and, in turn, change the nearshore hydrodynamics (e.g., Osorio-Cano et al., 2019;Reguero et al., 2019). Fiji, like many Pacific Island countries has reported coastal degradation over recent years. A study investigating the benthic habitat of the reef at a nearby site on the Coral Coast estimated about ∼48% of algae cover, ∼19% of sand and gravel cover, and only ∼7.5% of coral cover on the reef flat (Singh et al., 2019). It is, therefore, feasible that a developed coastal area, such as the Coral Coast, would have reduced friction coefficients compared to more pristine and less degraded coral reefs (e.g., Van Dongeren et al., 2013;Rogers et al., 2017).

Empirical Models
Because many large-scale studies of coastal flooding and extreme TWLs employ empirical formulations to derive η wave and maximum TWLs, we compared the multiple linear regression model and the XBeach simulations with two empirical solutions. We applied the same methodologies as Vitousek et al. (2017) and Vousdoukas et al. (2018) to reproduce the η wave and TWL.
while Vitousek et al. (2017) calculate η wave using the empirical relationship for the 2% exceedance run-up on dissipative beaches (Stockdon et al., 2006): where L 0 is the wavelength obtained using linear wave theory for deep water L 0 = gT 2 /(2π).
Both define the maximum TWL as

Offshore Waves
To understand the nature of the May 2018 swell, the atmospheric and wave conditions in the swell-generation zone were analyzed.  Figure S2). The strong winds were the result of a low-pressure system (<980 hPa minimum pressure) tracking south of New Zealand while a large stable high-pressure system (>1030 hPa maximum pressure) covered the extent of Australia. Maximum H s in the swell-generation area exceeded 6 m while peak periods of the northward propagating swell exceeded 15 s (Supplementary Figure S3). The swell reached the Coral Coast around 22:00 (UTC) on May 24, 2018, when the wave buoy recorded an increase in H s , T m , and T P . The swell was characterized by two distinct peaks in H s . The first peak in H s (4.67 m) occurred around 9:45 (UTC) on May 25, and the second peak (4.91 m) was recorded at 11:20 (UTC) on May 27 (Figure 2A). As such, the maximum H s exceeded the 99th percentile H s (Bosserelle et al., 2015b) by 25%. The maximum observed T P exceeded 17 s. The swell direction remained constant from the southwest between 190 • and 220 • throughout the entire event. The maximum recorded single wave height (derived from the vertical water displacement) was 9.32 m at 10:21 (UTC) on May 25.

Nearshore Waves and Water Levels
The largest instantaneous TWL was recorded during the second peak of the swell on May 27, 2018, at high tide (2.03 m above MSL, Figure 2B). Consistent with previous findings (e.g., Becker et al., 2016) IG and short wave heights followed the tidal signal with the largest wave heights during high tides and smallest wave heights during low tides (Figures 2C,D). On the reef flat, the height of IG waves exceeded the height of short waves significantly throughout the entire period with the maximum IG H s of 0.79 m during the second peak of the swell. The maximum H s of the short waves was 0.35 m ( Figure 2E). η wave also followed a tidal signal with the maximum setup observed during low tides and lower η wave during high tides. The maximum η wave was 0.91 m and occurred during the second peak of the swell. The out-of-phase signals of η wave and wave heights resulted in similar TWL increases during high and low tides. The max 2% total non-tidal water level was 1.42 m, and the maximum observed non-tidal water level was 1.64 m above η tide and occurred during low tide. During high tide, when the η wave decreased to 0.25 m, the observed non-tidal water level was still 1.62 m above η tide due to the large IG and short waves (Figure 2F).
The contribution of the different components (tide, short waves, IG waves, and setup) to the maximum instantaneous TWL was examined by calculating the percentage of each component to the hourly maximum observed water level ( Figure 2G). Between May 24 and 29, 2018, tides on average contributed to 55% of the maximum recorded TWL, η wave to 25%, IG waves to 11%, and short waves to 9%; however, the contribution varied significantly. Unsurprisingly, before the swell reached Fiji (before May 25, 2018), the TWL was largely tidally dominated (85-90%) as the offshore waves were small. As H s,inc increased with the arrival of the swell, the contribution of the tides to the TWL decreased and the contribution of setup, short, and IG waves increased. Particularly during low tides, η wave contributed the most to the maximum TWL while short and IG waves superseded the contribution of η wave during high tide. On May 27, even at high tide, η tide only contributed to ∼47% of the maximum TWL while η wave (10%), IG waves (24%), and short waves (19%) accounted for the rest.  To better explain these relationships, we performed a multiple linear regression analysis where offshore H s and η tide were the independent variables and where H s,SW , H s,IG , η wave , and hourly maximum instantaneous TWL were the dependent response variables (Figure 3). We found highly statistically significant (pvalues < 0.05) linear relationships between the predictor and the response variables: H s,SW R 2 = 0.88 , H s,IG R 2 = 0.86 , η wave R 2 = 0.92 , and the hourly maximum instantaneous TWL R 2 = 0.91 . This implies that 91% of the hourly maximum TWL near the shore can be explained by the incident wave heights and the tides.

Numerical Model Results
The XBeach model was able to accurately reproduce the nearshore water levels and waves (RMSE = 0.12 m) (Figure 2B). H s,IG (Figure 2C) and the η wave ( Figure 2D) were captured well with RMSEs of 0.07 and 0.12 m, respectively; however, the peaks in η wave were slightly underestimated. Simulated H s,SW were systematically smaller than the observed (RMSE of 0.16 m). This underestimation in η wave and H s,SW is likely responsible for the overall underestimation of TWLs (Figure 2B).
We performed a multivariate sensitivity analysis (Supplementary Figures S4, S5) to evaluate the model response to changes in friction coefficients. Simulated η wave was mainly sensitive to the wave friction coefficient while IG wave heights were sensitive to both bottom friction coefficient for flow and short wave friction coefficient (Supplementary Figures S4, S5). In comparison, wave heights and η wave were less sensitive to β and γ . The simulations were used to assess the inundation extent and showed that the flooding was mainly limited to the immediate coastline and only few inhabited areas (Supplementary Figure S6A). This estimation agrees with verbal witness accounts that pointed out water depth and water marks on land and stated most damage and impact during the swell occurred on buildings in immediate proximity to the ocean. The simulated inundation extent was calculated by summing up all inundated grid points. Throughout the swell, an area of 67,620 m 2 within the model domain was flooded.
The simulation with increased mean sea levels suggests the inundated area would be ∼97% larger (133,110 m 2 ) and include large areas of the coastal settlements under the RCP 8.5 scenario (Supplementary Figure S6B). Although the simulation did not account for future changes in wave climate, coral cover, bathymetry, topography, or potential adaptation measures, it serves as a first estimate of changing hazards due to climate change. The maximum H s,IG increased by 28% while the maximum η decreased by ∼61%. H s,SW also increased (by 244%); however, as the model failed to accurately simulate H s,SW under current sea levels, the projected changes might be inaccurate.

Comparison of Different Approaches to Estimate TWL
Comparison between the empirical formulas, the multiple linear regression model, and the XBeach simulations (Figure 4) show clear limitations of the empirical formulations. Setup was tidally modulated, and the parametric equations, which are functions of just H s,inc , were not able to capture the tide-induced variability in setup. This resulted in a significant overestimation of η wave during low tides as well as during most high tides ( Figure 4A). Although the TWL max predictions were better than the predictions of η wave , maximum water levels were too low compared to the observations ( Figure 4B). The XBeach simulation performed better than the empirical formulations and was able to capture η wave and maximum water levels well despite slightly underestimating the TWL (Supplementary Table S7).
The multiple linear regression model was able to reproduce the observations most closely, which is not surprising as the model was trained with the observations from this particular event.

DISCUSSION AND CONCLUSION
In the present study, we investigated a large swell in May 2018 that was generated by an extra-tropical low-pressure system in the Southern Ocean and Tasman Sea. The waves generated by the system propagated through the South Pacific and caused extreme TWLs and coastal inundation along Fiji's Coral Coast. The maximum TWL near the shore was caused by several compounding mechanisms, namely η tide (55% average contribution to the maximum TWL), η wave (25% average contribution to the maximum TWL), IG waves (11% average contribution to the maximum TWL), and short waves (9% average contribution to the maximum TWL). Consistent with previous findings (e.g., Péquignet et al., 2011;Pomeroy et al., 2012;Becker et al., 2016Becker et al., , 2014 results indicated that H s,SW , H s,IG , and η wave on the reef were modulated by H s,inc and η tide . Consequently, the maximum hourly instantaneous TWL could be described as a function of η tide and H s,inc . Simulations of the swell with the numerical model XBeach were able to reproduce the main mechanisms contributing to the extreme TWL and the inundation extent. Simulations of the same event under the RCP8.5 sea-level scenario suggest an increase in inundation area by 97%. Although the contribution of η wave to the maximum TWL decreased (by 61%), the contribution of waves (both short and IG waves) increased.
A comparison between a multiple linear regression model, XBeach, and two empirical equations to estimate η wave and maximum TWL Vousdoukas et al., 2018), reveals that both XBeach and the linear regression model were better suited to reproduce the nearshore water levels than the empirical formulations.
The benefit of the linear regression model is its ease of use as it does not rely on computationally expensive simulations. The linear relationship between H s,inc , η tide , and nearshore TWL could be a useful tool to efficiently forecast extreme TWLs along the Coral Coast. However, it does not provide any information about inundation or wave run-up. Furthermore, although the model performed well for the event with which it was trained, it is yet to be tested against different swells (e.g., different wave period or direction).
The numerical model XBeach, on the other hand, is well established and tested in fringing reef environments (Van Dongeren et al., 2013;Buckley et al., 2014). Although XBeach is a useful tool to examine single events, it is less suitable as an operational forecast system due to the high computational cost. However, recent advances in hybrid-downscaling of coastal flood models (e.g., Camus et al., 2011;Rueda et al., 2019Rueda et al., , 2016 could provide an adequate solution for operational forecast systems in the region. In this study, the model slightly underpredicted H s,SW and the maximum TWL. Although the best available bathymetry and topography data were used to set up the model, some of the model error may still be attributed to the lack of better baseline data (e.g., Lidar topography and bathymetry). Another possible solution to improve the model could be to apply a varying friction coefficient depending on marine habitat across the model domain (Passeri et al., 2018). Coral reefs are complex systems in terms of geometry and bottom roughness; hence, incorporating this complexity into the numerical model might improve the results (e.g., Osorio-Cano et al., 2019;Reguero et al., 2019;Silva et al., 2019).
The results presented here highlight the need for customized, locality-specific coastal hazard assessments and inundation forecast systems based on observational data in the South Pacific as the commonly applied empirical equations often fail to capture the complex physical processes within reef environments. Localized inundation thresholds (Smith and Juria, 2019), downscaled numerical models, and in situ wave and sealevel observations are critical tools to understand swell-driven hazard in fringing reef environments. Decision makers are urged to invest in regional coastal monitoring programs, including long-term observations and downscaled models, to find adequate solutions to mitigate the inundation risks Pacific Island nations will face over the coming decades.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
MW conceived the study and performed most of the simulations, analysis, and writing and some of the field data collection. JA assisted with field work and data analysis. AE assisted with writing and data analysis. NJ assisted with field work and provided the model grid. AD assisted with field work and writing. HD contributed data and other information and assisted with modeling, writing, and analysis. All authors contributed to the article and approved the submitted version. FUNDING MW, AE, NJ, and HD were funded through the Coastal Inundation Forecast Demonstration Project (CIFDP) funded by the World Meteorological Organization (WMO) and the Korean Government and through the Climate and Oceans Support Program in the Pacific (COSPPac) funded by the Australian Government. AD was funded through a JEAI grant from the French Research Institute for Sustainable Development (IRD) to the University of the South Pacific.