Global mangrove distribution on the open coast is controlled by waves

Site selection

This study is based on 20 years of data (3-h temporal resolution) from 400 globally distributed sites. Sites were selected from all continents containing mangroves (Africa, Asia, North America, Oceania and South America) and chosen to extend to their latitudinal boundaries, therefore encompassing a broad range of genera and climates. The selected sites included 200 with mangroves on the open coast and 200 sites without mangroves (Fig. 2). Sites with and without mangroves were identified on the basis of distribution data from Global Mangrove Watch⁵⁰ and were cross-referenced with Google Earth satellite imagery from 2025. Sites with mangroves present were selected only where mangroves seemed to be exposed to waves, based on visual inspection of the satellite imagery. As such, back-barrier mangroves protected from offshore wave conditions due to the presence of a dune or reef were not included in the dataset. Sites protected by a barrier reef or an island at an offshore distance greater than the resolution of the mesh of the wave model (about 20 km) were permitted in the analysis (for example, sites near the Great Barrier Reef). Embayed sites were permitted where the body of water (for example, bay or gulf) had a width greater than 20 km.

To reduce the influence of climate factors on mangrove distribution, non-mangrove sites were selected in marine ecoregions that contain mangroves. These non-mangrove sites were identified using 2025 Google Earth satellite imagery in coastal environments that seem to comprise sand or mud (rocky headlands, for example, were excluded). Sites with any visible development in the satellite imagery that could preclude mangrove growth such as hardened shorelines and seawalls were not selected for the analysis.

Data collection

This study was conducted using wave data from the Global Ocean Waves Reanalysis⁴⁸ (https://doi.org/10.48670/moi-00022) of the EU Copernicus Marine Service Information. This dataset contains spectral wave data at 0.2° × 0.2° spatial resolution globally and 3-h temporal resolution. The product is based on MFWAM, a third-generation spectral wave model. This model calculates the wave spectrum using ERA5 climatological winds⁵¹, GLORYS12V1 currents⁵² and ETOPO2 bathymetry⁵³. Twenty years of wave data from 1 January 2003 to 1 January 2023 were used in the analysis. Three swell partitions were considered: the primary, secondary and wind components of the spectrum. The bearing of the coastline at each site was measured to the nearest 5° in QGIS with Google Earth basemap (map data 2025 Google).

WOOs

Significant wave height and direction were chosen to represent the magnitude of wave interaction at the shoreline. The physical mechanisms leading to mangrove mortality as a result of waves are forces due to near-bed orbital flow velocity and erosion due to bed shear stress¹⁹. For mangrove seedlings (which have diameters D ≪ 1 m), inertial forces due to near-bed flow acceleration can be considered negligible. The drag force F_D due to flow velocity (u) acting on a submerged structure is given by

$${F}_{{\rm{D}}}=\frac{1}{2}{C}_{{\rm{D}}}\,\rho {Au}|u|$$

(1)

where C_D is a drag coefficient, A is the projected area and ρ is seawater density in this case.

Mean bed shear stress is given by

$${\tau }_{{\rm{b}}}=\frac{1}{2}{f}_{{\rm{w}}}\rho {u}^{2}$$

(2)

where f_w is the wave friction factor. Both physical mechanisms for mangrove seedling dislodgement—drag forces and erosion—are therefore driven largely by the bed flow velocity.

The maximum horizontal component of the bed flow velocity due to waves is given by

$$u=\frac{\omega H}{2\sinh (kh)}$$

(3)

where H is the wave height, ω is the wave angular frequency ($\omega \equiv \frac{2\pi }{T}$), T is the wave period, h is the water depth and k is the wavenumber ($k\equiv \frac{2\pi }{\lambda }$). In shallow water, where wavelengths (λ) are approximated by $\lambda \approx T\sqrt{{gh}}$ and $\sinh (kh)\approx kh$, this simplifies to

$$u=\frac{H\sqrt{g}}{2\sqrt{h}}.$$

(4)

where g is gravitational acceleration. This approximation is used to justify the choice of the significant wave height H_s as the offshore wave property most appropriate to represent mangrove seedling establishment thresholds.

As an offshore wave approaches the coast, it is transformed by refraction. The refracted wave height H_r is given by

$${H}_{{\rm{r}}}={K}_{{\rm{r}}}{H}_{0}$$

(5)

where H₀ is the incident wave height and K_r is the refraction coefficient given by

$${K}_{{\rm{r}}}=\sqrt{\frac{\cos {\theta }_{0,\perp }}{{\cos \theta }_{\perp }}}$$

(6)

and ${\theta }_{0,\perp }$ is the angle relative to shore normal of the offshore wave and ${\theta }_{\perp }$ is that of the refracted wave.

It was assumed that waves at the shoreline travel perpendicular to the shore, so that the wave height near the shore can be expressed as

$${H}_{{\rm{r}}}={H}_{0}\sqrt{\cos {\theta }_{0,\perp }}.$$

(7)

To capture the component of the wave spectrum most impactful to the coastline for each 3-h timestep, the EWH at time t was therefore defined as

$${\rm{E}}{\rm{W}}{\rm{H}}(t)=\mathop{\text{max}}\limits_{i\in \text{Primary},\text{Secondary},\text{Wind}}\{{{H}_{{\rm{s}},}}_{i}(t)\cdot \sqrt{\cos {\theta }_{\perp ,i}(t)}\}$$

(8)

where H_s is the offshore significant wave height and ${\theta }_{\perp }$ is the mean wave direction relative to the beach normal for each of the spectrum partitions considered. This represents the component of the wave spectrum that is expected to have the most impact at the site at time t. The maximum of the swell partitions was used, rather than bulk wave statistics of the whole spectrum, to better capture scenarios in which the spectral mean wave direction would result in a low value of EWH (${\theta }_{\perp ,m}\approx \pm {90}^{^\circ }$), but a particular swell partition is approaching at an angle that is more impactful at the shore (${\theta }_{m}\approx {0}^{^\circ }$). A diagram demonstrating the calculation of this parameter is shown in Extended Data Fig. 1. This defined parameter accounts for the incidence of non-locally generated swell and is temporal, enabling the identification of low-energy windows key to mangrove establishment.

Twenty years of spectral wave data were extracted from the nearest offshore non-null node of the wave dataset⁴⁸ for each of the 400 analysis sites. Wave nodes were classified as offshore if the bearing between the node and the analysis site approached the shore at an angle between −90° and 90° relative to shore normal. In doing so, three mangrove and three non-mangrove sites were removed from the datasets as a result of the surrounding offshore nodes all being null. At each site for each 3-h timestep, the EWH was calculated.

A deterministic model to predict mangrove presence was then created using only the training sites. The establishment thresholds for mangrove seedlings were calculated on the basis of wave data from the sites with mangroves, using the following process. Time periods representing mangrove seedling establishment up to 1 year were defined as D_s = {1 day, 2 days, 5 days, 10 days, 2 weeks, 4 weeks, 8 weeks, 12 weeks, …, 52 weeks}. These increasing durations were selected to represent the increasing tolerance of mangrove seedlings to disturbance as they grow³². This set was defined with higher resolution for smaller durations to capture the rapid initial anchoring of seedlings²² and durations up to 1 year were considered to capture the seasonality of sites. At each site, the minimum thresholds that have a period of non-exceedance for the first critical duration of 1 day (1d)—were determined for each month (M) in the dataset. These thresholds C_M,1d were defined as

$${C}_{M,1d}=\text{min}\{c:{\rm{\exists }}\,{t}_{0}\in M:{{\rm{E}}{\rm{W}}{\rm{H}}}_{M}(t) < c\,{\rm{\forall }}\text{t}\in [{t}_{0},{t}_{0}+1\,\text{day}]\}$$

(9)

where EWH_M(t) includes all values of the EWH for the given month at the site. Each of these durations was then considered a possible candidate for the early establishment of a mangrove propagule. The maximum value of the EWH for co-originating durations $D\in {D}_{{\rm{s}}}$ was then determined. These durations represent the windows in which hydrodynamic thresholds for a seedling of that age should not be exceeded, for successful establishment. The non-exceedance threshold for each duration D following the 1-day threshold C_M,1d for each month M was defined as C_M,D and calculated as

$${C}_{M,D}=\text{max}\{{\rm{EWH}}(t)\}\text{for t}\in [{t}_{0},{t}_{0}+D]$$

(10)

As the initial thresholds C_M,1d were calculated from the beginning of each month, this resulted in 12 sets of non-exceedance values (each containing one threshold per duration $D\in {D}_{{\rm{s}}}$) per year at each site. Each of these sets was then considered a candidate for a seedling establishment event. In reality, propagule availability is a bottleneck to establishment in many environments¹⁵. However, global trends have shown that propagule abundance is typically minimal during winter⁵⁴ when time waves are more energetic. This suggests that, globally, optimal establishment events in terms of wave forcing will typically occur when propagules are available. There are, however, exceptions to this trend and future site-specific studies should therefore take local propagule seasonality into account.

The sites allocated for training the model (all sites excluding the one on which the model is being assessed, as per the leave-one-out cross-validation (LOOCV) method; see ‘Model performance’) were used to determine the establishment thresholds of mangroves. Sites with and without mangroves had different distributions of non-exceedance values (Fig. 3a). A series of Anderson–Darling k-sample tests (see ‘Model performance’) were conducted to compare these distributions for each duration, with the results shown in Extended Data Table 1. To ensure data independence, the average non-exceedance value for each duration at each site was used (n₁ = n₂ = 197). The tests indicate that the null hypothesis that mangrove and non-mangrove sites are drawn from the same distribution can be rejected for all durations at the 2% significance level after Bonferroni correction for 18 tests (as P < 0.0011 for all tests and α_corrected = 0.11% is the per-test significance level corresponding to a family-wise error rate of 2%; see ‘Model performance’).

The 65th percentile of the threshold values at sites with mangroves for each duration were selected as the upper limits of establishment for a seedling of that age. This percentile was chosen by determining the value that optimized (≈1) the ratio of exceedance at the non-mangrove sites and non-exceedance at the mangrove sites (that is, the values that resulted in the best performance for specificity and sensitivity). This resulted in a hydrodynamic habitat threshold C_D corresponding to each duration $D\in {D}_{{\rm{s}}}$ (Fig. 3b). These thresholds were then recalculated using only sites in micro-, meso- and macro-tidal environments to determine how tidal range may affect hydrodynamic tolerance (tidal classification from a dataset³⁰ based on the FES2014 global tide model³¹).

The alignment of these calculated thresholds with previous empirical studies^27,32 was then assessed. The present work suggests a threshold of 1.13 m for the EWH for 12-month-old seedlings (Fig. 3; see ‘WOOs’, Methods for calculation). If the shallow water depth breaking limit is taken to be $\frac{H}{h}=0.78$ (ref. ⁵⁵) (where H is wave height and h is water depth), this wave can be assumed to be broken when the depth becomes shallower than 1.44 m, beyond which linear wave equations will no longer apply. In shallow water, the resulting near-bed wave orbital velocity will be u = 1.46 m s⁻¹ (equation (4)). It has previously been estimated³² that seedlings (less than 1 year old) may become susceptible to uprooting at flow velocities 1.2–1.5 ms⁻¹; this aligns with the present result.

If a standard friction factor of f_w = 0.01 (ref. ⁵⁶) and saltwater density of ρ = 1,025 kg m⁻³ are used, the bed shear stress corresponding to this EWH limit is τ_b = 10.88 Nm⁻² (equation (2)). This lies within 2 s.d. of the shear-stress limit that has been imposed previously²⁷ in models for seedlings of the same age.

The following predictive approach was then applied to all sites (197 with and 197 without mangroves). It was asserted that, for a mangrove seedling to establish in a particular year, all hydrodynamic habitat thresholds should be satisfied. Thus, if the following condition is satisfied at a given site for a particular year, that year is considered suitable for mangrove establishment.

$${\rm{\exists }}\,{t}_{0}:{\rm{EWH}}(t) < {C}_{{\rm{D}}}\,{\rm{\forall }}t\in [{t}_{0},\,{t}_{0}+D]{\rm{\forall }}\,D\in {D}_{{\rm{s}}}\,$$

(11)

Previous studies^17,18 have indicated that a single set of WOOs is sufficient to initiate mangrove establishment, meaning that mangrove stands may be present at sites despite experiencing unsuitable conditions outside the WOO. It was therefore asserted that, if the previous condition was met for any of the years considered, the site could be considered suitable for mangroves. The resulting model output is a binary result predicting the presence or absence of mangroves at each site. The model results were first assessed for all sites collectively, then for sites with Avicennia, Rhizophora or Sonneratia reported on the seaward fringe³³ to determine how different mangrove genera may respond to hydrodynamic forcing.

Model performance

The performance of the predictive model developed in this work is assessed using accuracy, sensitivity and specificity. These quantities were calculated as described in the following expressions.

$${\rm{Accuracy}}\equiv \frac{\text{Number of true positives}+\text{number of true negatives}}{\text{Total number of points}}$$

(12)

$${\rm{Sensitivity}}\equiv \frac{\text{Number of true positives}}{\text{Number of true positives}+\text{number of false negatives}}$$

(13)

$${\rm{Specificity}}\equiv \frac{\text{Number of true negatives}}{\text{Number of true negatives}+\text{number of false positives}}$$

(14)

True positives are defined as points that are predicted to have mangroves and do, whereas false positives are defined as points that are predicted to have mangroves but do not. True negatives are those predicted to be without mangroves that do not, whereas false negatives are those predicted to be without mangroves that do have mangroves.

The model performance was validated using the LOOCV method⁵⁷. LOOCV involves assessing the model prediction at all sites in the study. For each point, all other points are used to train the model, and the trained model is then used to predict the outcome (mangrove or non-mangrove) for the remaining point. In this case, the training of the model means that the 65th percentile of the non-exceedance value of EWH at all other sites was used as the establishment threshold to predict mangrove establishment at the remaining site.

The statistical significance of the difference between distributions of non-exceedance values of the EWH at sites with and without mangroves was assessed using the Anderson–Darling k-sample test and corrected for several tests using the Bonferroni correction. The Anderson–Darling k-sample test⁵⁸ is a modification of the Anderson–Darling test⁵⁹ and is a non-parametric statistical test to evaluate the null hypothesis that two or more groups are drawn from the same underlying distribution. This test was selected because it does not require specification of the underlying population distribution and it is sensitive to differences across the entire distribution. To ensure data independence, the mean non-exceedance value (calculated across all sampled monthly starting points) for each duration at each site was used. This test computes the test statistic (A²), which can be compared with a critical value to determine whether the null hypothesis can be rejected at the corresponding significance level. When conducting several statistical tests, the probability of a type-I error (incorrect rejection of a true null hypothesis) occurring in at least one of the tests is greater than the probability of this error occurring in any one test individually. To adjust for this, a Bonferroni correction was applied, which defines the adjusted significance level as⁶⁰

$${\alpha }_{{\rm{corrected}}}=\frac{\alpha }{m}$$

(15)

where m is the number of tests. In this study, tests were considered statistically significant at the Bonferroni-adjusted significance level corresponding to a family-wise error rate of 2%. This corresponds to a significance level for each test of α_corrected= 0.11%.

Case studies

Two case study sites were chosen that had a coastline with both mangrove and non-mangrove regions. The first was selected in the Palk Strait of the Indian Ocean (9.599069° N, 78.911753° E; 9.238059° N, 79.286306° E); the second in the Gulf of Thailand of the Pacific Ocean (9.342269° N, 99.187294° E; 9.221632° N, 99.352620° E). Equidistant points (0.05° and 0.1° spacing) were created along the coastline selected for the case study. The shoreline at each of the points was then visually classified as mangrove, non-mangrove or developed. Mangrove sites were classified as described previously (‘Site selection’, Methods) and points were classed as developed where houses or structures were built along the shoreline. The model was then applied as described previously (‘WOOs’, Methods), using thresholds calculated using all valid mangrove sites in the global analysis (197). Each of the points was then classed as mangrove suitable or unsuitable. This was repeated using the historical values of EWH decreased and increased across a range of percentages.

Global sensitivity to wave climate

The sites that were correctly identified by the model (as described in ‘WOOs’, Methods) as mangrove or non-mangrove were reassessed under altered conditions to assess their sensitivity to changes in wave climate. The 20 years of historical data used to predict current presence were multiplied by factors from 55% to 145% (changes of −45% to +45%). The suitability of the sites for mangroves was then reassessed using these data.