Coral reef structural complexity loss exposes coastlines to waves
Ecological sampling and structural complexity profilesThe ecological sampling consists of 10 surveys, taking place in 2005 and from 2008 to 2016, and documents changes in coral colony abundance and size distributions (i.e. width, length, and height) for the three most conspicuous taxa (i.e. Acropora, Pocillopora, and Porites) within a 10 m2 transect on the outer slope23. To quantify reef structural complexity, we built a 3D model of the coral assemblages distributed along a cross-section of the reef substrate separating the 20 m water depth from the reef crest, representing a 160 m stretch along the reef slope (Fig. 1). First, we take 200 overlapping high-resolution photos (300 dpi) of 10 individual corals from each species (i.e. n = 30 coral colonies) and built 3D models using the Agisoft Metashape software24, capturing intra- and inter-species morphological variability (Fig. 1). Then, we systematically and randomly select one of the ten 3D coral models for each taxon to add to the substrate until that the sum of the planar area for each 3D coral models match with the coral cover reported for each taxon and for each year23. We randomly place coral colonies along the 160 m reef cross-section going from 20 m depth to the reef crest (Fig. 1). The individual coral 3D models are resized in width, length, and height according to ecological surveys, and, randomly rotated between − π/2 and π/2 to ensure ecological variability. Finally, we estimated structural complexity of the 3D coral assemblage model using the function rumple_index of the LidR package25 in R 4.0.026. We repeat this approach 100 times for each year, resulting in a total of 1000 reef structural complexity profiles. Our estimates are consistent with previous reef structural complexity estimates at this location27.Figure 1(a) Representation of the three different coral species (Acropora hyacinthus in red, Pocillopora cf. verrucosa in yellow, and Porites lutea in blue). (b) A representaitive Ha’apiti reef cross-section simulation (one of 1000 total simulations) on the outer slope across a water depth range of 0–20 m.Full size imageHydrodynamic and topographic measurementsMo’orea (French Polynesia) is encircled by coral reefs, 500–700 m wide with a dominant swell direction coming from the southwest. In this study, we focus on Ha’apiti, a site with a southwest orientation that is considered as a high-energy site28. We extract 30-year offshore wave data (1980–2010) from a wave hindcast8,29 (Fig. 2a). We also collect high-frequency, in situ wave data using INW PT2X Aquistar and DHI SensorONE pressure transducers (PTs), which are logged at 4 Hz30. The sensors are installed at four locations along a cross-shelf gradient (Fig. 2b,c) covering a 250 m long stretch, including sections through the fore reef, reef crest, and reef flat. Pressure records are corrected for pressure attenuation with depth31 and are split into 15-min bursts30.Figure 2(a) Histogram of the offshore wave height (m) at Ha’apiti, Mo’orea (French Polynesia) in 2016. (b) Aerial view of Ha’apiti (WorldView-3 imagery) with an outline of the wave transect and sensor location. The ecological sampling took place near the S1 location c. Topographic cross-section of the wave transect and position of the sensors on the sea bottom.Full size imageThe beach profile and the reef morphology are measured using airborne bathymetric and topo-bathymetric lidar surveys conducted in June 2015 by the Service Hydrographique & Océanographique de la Marine (SHOM). The bathymetric data are defined by the combination of bathymetric laser (for the submerged part of the beach) and topo-bathymetric laser (for the subaerial beach). The data come at 1 m resolution and are available at https://diffusion.shom.fr.Hydraulic roughness vs structural complexitySpectral attenuation analysis of the water level measurements32,33 is used to estimate the Nikuradse (hydraulic; kn) roughness34 of the coral reef surface along the beach profile sections covered by the pressure transducers. The method is described in detail in the references provided above and uses the conservation of energy equations to obtain estimates of wave energy dissipation from friction. We obtain more than 300 kn estimates for each pair of sensors, each representing a different geomorphologic section. Since the field measurements took place in 2015, the kn outputs obtained from the fore reef section concur with the reef structural complexity estimates of that year (Fig. 3). Then, we define a coefficient factor according to the geomorphologic section as ⍺back reef = kn, back reef/kn, fore reef and ⍺reef crest = kn, reef crest/kn, fore reef. We carefully delineate the sandy section from the reef sections within the cross-shelf gradient (i.e. within the reef flat, lagoon section) and apply the following procedure. First, for the reef sections, we apply the relationship between the reef structural complexity and kn (Fig. 3) to convert our reef structural complexity estimates into continuous kn profiles through Monte Carlo simulations, using the coefficient factor of each geomorphologic section (e.g., forereef, reef crest, and back reef). Second, for the sandy section, we define kn on the grounds of the mean grain size (d50 = 63 μm). Applying this workflow (Fig. 3), we obtain 100 continuous kn profiles for each year (i.e. n = 1000 kn profiles in total).Figure 3Flow chart illustrating how the kn profiles have been obtained along the cross-section at Ha’apiti. The relationship between the Structural complexity (SC) and the Nikuradse roughness (kn) measurements can be described as kn = 0.01 × SC2.98.Full size imageHydrodynamic modelNearshore wave propagation is simulated using a nonlinear wave model based on the Boussinesq Equations35. The rationale of using a Boussinesq type model instead of other types of models (e.g. SWAN) is that the former is able to describe in detail (i.e. 1 m grid resolution) several hydrodynamic parameters (e.g. nearshore nonlinear wave propagation, shoaling, refraction, dissipation due to the bottom friction and breaking and run-up) in the swash zone. The model is defined as follows:$$frac{partial U}{partial t}+frac{1}{h}frac{partial {M}_{u}}{partial x}-frac{1}{h}Ufrac{partial left(Uhright)}{partial x}+gfrac{partialupzeta }{partial x}=frac{left({d}^{2}+2partialupzeta right)}{3}frac{{partial }^{3}U}{partial {x}^{2}partial t}+{d}_{x}hfrac{{partial }^{2}U}{partial xpartial t}+frac{{partial }^{2}}{3}left(Ufrac{{partial }^{3}U}{{partial x}^{3}}-frac{partial U}{partial x}frac{{partial }^{2}U}{partial {x}^{2}}right)+dfrac{partialupzeta }{partial mathrm{x}}frac{{partial }^{2}U}{partialupzeta partial mathrm{t}}+d{d}_{x}Ufrac{{partial }^{2}U}{partial {x}^{2}}+{d}_{x}frac{partialupzeta }{partial mathrm{x}}frac{partial mathrm{U}}{partial mathrm{t}}-dfrac{{partial }^{2}}{partial mathrm{x}partial mathrm{t}}left(delta frac{partial mathrm{U}}{partial mathrm{x}}right)+E-frac{{tau }_{b}}{rho h}+B{d}^{2}left(frac{{partial }^{3}U}{partial {x}^{2}}+gfrac{{partial }^{3}upzeta }{partial {x}^{3}}+frac{{partial }^{2}left(Ufrac{partial U}{partial x}right)}{partial {x}^{2}}right)+2Bd{d}_{x}left(frac{{partial }^{2}U}{partial xpartial t}+gfrac{{partial }^{2}upzeta }{partial {mathrm{x}}^{2}}right),$$
(1)
where, U is the mean over the depth horizontal velocity, ζ is the surface elevation, d is the water depth, uo is the near bottom velocity, h = d + ζ, ({M}_{u}=left(d+zeta right){u}_{0}^{2}+delta ({c}^{2}-{u}_{0}^{2})), δ is the roller thickness determined geometrically36, E is an eddy viscosity, τb is the bed friction term and B = 1/1535.In this work the wave breaking mechanism is based on the surface roller concept36. However, in the swash zone, surface roller is not present and the eddy viscosity concept is used to describe the breaking process. The term E in Eq. (1) is written:$${mathrm{E}}_{{mathrm{b}}_{mathrm{x}}}= {mathrm{B}}_{mathrm{b}}frac{1}{mathrm{h}+upeta }{left{{{mathrm{v}}_{e}left[left(mathrm{h}+upeta right)mathrm{U}right]}_{mathrm{x}}right}}_{mathrm{x}},$$
(2)
where ({v}_{e}) is the eddy viscosity coefficient:$${mathrm{v}}_{mathrm{e}}={{ell}}^{2}left|frac{partial {mathrm{U}}}{partial {mathrm{x}}}right|,$$
(3)
where ({ell}) is the mixing length ({ell}) = 3.5 h και Βb37.The width of the swash zone is assumed to extend from the run-down point (seaward boundary) up to the run-up point (landward boundary). We start from a first estimate of the run-up R using the Stockdon formula38 and the depths below R/4 are considered as the swash zone, using Eq. (2). The final wave run-up height R which comes as output is estimated by the model.The ‘dry bed’ boundary condition is used to simulate run-up35. The numerical solution is based on the fourth-order time predictor–corrector scheme39. Therefore, the bed friction term τb is calculated such as:$${tau }_{bx}=frac{1}{2}rho {f}_{w}Uleft|Uright|,$$
(4)
where fw is the bottom friction coefficient40, which is an explicit approximation to the implicit, semi-empirical formula given by Jonsson, 196741.$${f}_{mathrm{w}}=mathrm{exp}left[{5.213left(frac{{mathrm{k}}_{mathrm{n}}}{{mathrm{alpha }}_{0}}right)}^{0.194}-5.977right],$$
(5)
where αo is the amplitude of the near-bed wave orbital motion and kn is the Nikuradse roughness height.Simulations and post processingWe use our wave propagation model to assess how different coral reef states affect the impact waves have on the coast. We run an ensemble of 10,000 simulations that covers all the possible combinations of (i) 10 bottom roughness profiles expressing the different observed coral reef states (i.e. healthy vs. not unhealthy); and (ii) 1000 percentiles of wave conditions. The wave conditions are produced as follows: (i) from the weekly values, we estimate all significant wave height (Hs) percentiles from 0.1 to 100, with a step of 0.1; (ii) the resulting 1000 Hs values are linked to the corresponding peak wave period Tp using a copula expressing the dependence of the two variables42. The output of the simulations is the nearshore Hs and 2% exceedance run-up (R2%) height for each of the 1000 conditions and 10 coral reef states. To quantify how the coral reef states are altering wave propagation during extreme events, we apply extreme value analysis to estimate the R2% for different return periods43. We then compare how the return period curves changed from the two coral reef states and we define the change in frequency of extreme R2% under unhealthy coral reefs. It is important to highlight that the tidal range is More
