Coral reef structural complexity loss exposes coastlines to waves

Ecological sampling and structural complexity profiles

The 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 m² transect on the outer slope²³. 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 software²⁴, 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 year²³. 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 package²⁵ in R 4.0.0²⁶. 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 location²⁷.

Hydrodynamic and topographic measurements

Mo’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 site²⁸. We extract 30-year offshore wave data (1980–2010) from a wave hindcast^8,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 Hz³⁰. 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 depth³¹ and are split into 15-min bursts³⁰.

The 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 complexity

Spectral attenuation analysis of the water level measurements^32,33 is used to estimate the Nikuradse (hydraulic; k_n) roughness³⁴ 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 k_n estimates for each pair of sensors, each representing a different geomorphologic section. Since the field measurements took place in 2015, the k_n 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} = k_{n, back reef}/k_{n, fore reef} and ⍺_{reef crest} = k_{n, reef crest}/k_{n, 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 k_n (Fig. 3) to convert our reef structural complexity estimates into continuous k_n 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 k_n on the grounds of the mean grain size (d₅₀ = 63 μm). Applying this workflow (Fig. 3), we obtain 100 continuous k_n profiles for each year (i.e. n = 1000 k_n profiles in total).

Hydrodynamic model

Nearshore wave propagation is simulated using a nonlinear wave model based on the Boussinesq Equations³⁵. 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:

where, U is the mean over the depth horizontal velocity, ζ is the surface elevation, d is the water depth, u_o 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 geometrically³⁶, E is an eddy viscosity, τ_b is the bed friction term and B = 1/15³⁵.

In this work the wave breaking mechanism is based on the surface roller concept³⁶. 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:

where ({v}_{e}) is the eddy viscosity coefficient:

where ({ell}) is the mixing length ({ell}) = 3.5 h και Β_b³⁷.

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 formula³⁸ 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-up³⁵. The numerical solution is based on the fourth-order time predictor–corrector scheme³⁹. Therefore, the bed friction term τ_b is calculated such as:

where f_w is the bottom friction coefficient⁴⁰, which is an explicit approximation to the implicit, semi-empirical formula given by Jonsson, 1967⁴¹.

where α_o is the amplitude of the near-bed wave orbital motion and k_n is the Nikuradse roughness height.

Simulations and post processing

We 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 (H_s) percentiles from 0.1 to 100, with a step of 0.1; (ii) the resulting 1000 H_s values are linked to the corresponding peak wave period T_p using a copula expressing the dependence of the two variables⁴². The output of the simulations is the nearshore H_s and 2% exceedance run-up (R_2%) 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 R_2% for different return periods⁴³. We then compare how the return period curves changed from the two coral reef states and we define the change in frequency of extreme R_2% under unhealthy coral reefs. It is important to highlight that the tidal range is < 20 cm at our study site, so will have negligible effect compared to other parameters. In systems with higher tidal range, tidal water level fluctuations should be also considered in the probabilistic framework.

Using the outputs of our simulations, we develop a Bayesian model which includes all the interdependencies between the run-up, the significant offshore wave height (H_{s, offshore}), and the reef structural complexity (SC). Our model is built in the R package brms^44,45 and follows the following structure:

where RU is the run-up (m), SC is the structural complexity, and H_s,offshore is the significant offshore wave height (m). The prior sampling is specified to follow a Gaussian ((mathcal{N}) (location, scale)) and a Gamma (Γ(shape, inverse scale)) distribution. We ran our models with four chains, 5000 draws per chain, and a warm-up period of 1000 steps, thus retaining 16,000 draws to construct posterior distributions. We verify chain convergence (n = 4) with trace plots and confirm that R_hat (the potential scale-reduction factor) is lower than 1.05⁴⁶. Drawing on our model, we predict the run-up according to six offshore wave height conditions (i.e. from 1 to 6 m with a step of 1 m) and the structural complexity gradient at Mo’orea from 2005 to 2016. In order to define the residual run-up height according to H_{s, offshore} and the reef structural complexity, we subtract from our run-up estimates the minimum value obtained for each offshore wave height condition.

Source: Ecology - nature.com