Case study: Kubulau District, Fiji

Villages and governments on Vanua Levu, the second largest island in the Fijian Archipelago, are becoming increasingly concerned that economic development (e.g., commercial agriculture, forest logging expansion) is impacting fisheries, tourism, and the ecological integrity of coral reefs and nearshore ecosystems61. For instance, a 20-year logging concession agreement was signed in 2006 for the Bua province on Vanua Levu, where Kubulau is located (Fig. 1). Growing demand for wood and root crops (i.e., taro and kava) for export increases pressure to convert native forest to pine forest in logging concessions, while commercial agriculture expands on land and slopes suitable for taro and kava monoculture. To this day, Kubulau District still has 70–80% of its native forest cover with relatively intact hydrologic connectivity between terrestrial, freshwater, and marine areas52,62, except for introduced pine tree plantations from previous logging activities along the coast52,62.

Figure 1 Study site. (a) Location of Kubulau District study site in Fiji. (b) Watershed boundaries, pour points, coral reef areas, and reef surveys are shown along with the customary fishing grounds (qoliqoli), periodically-harvested fisheries closures ( tabu ) and district level no-take MPAs. Full size image

To address these concerns, communities are currently working with government, NGOs and the private sector to design and implement ICM plans based on a national ICM frameworks for Fiji63. To address the tradeoffs inherent in designing ICM management interventions, information on priority areas for forest conservation or restoration that can promote coral reef ecosystem resilience under projected climate change is needed63. To support these efforts, we estimated the impact of projected land-use change on coral reefs under climate change scenarios in Kubulau District (population ~1,000) based on current and potential future logged areas. Our sediment model domain spanned all the watersheds discharging in traditional fishing grounds (qoliqoli), which include watersheds in Kubulau District and three large adjacent watersheds (Fig. 1b). In the customary fisheries management area (qoliqoli), coral reef habitats include fringing reefs, lagoons, patch reefs, and barrier reef extending over 30 km offshore. The local communities have implemented an MPA network of 21 community-managed periodically harvested fisheries closures (tabu areas) and three district-wide permanent no-take MPAs64.

Modeling approach overview

To spatially prioritize terrestrial conservation efforts across the landscape based on downstream coral reef impacts, we determined the impacts of projected climate change and land-use change on coral reefs and traced those back to the areas driving these impacts within each watershed. In order to do so, we coupled various land-use scenarios with an adapted linked land-sea modeling framework57 (Fig. 2). Previously developed for quantifying the effect of nutrient enriched groundwater on coral reefs, the modified framework links InVEST SDR version 3.259,60 to coral reef ecosystem state models calibrated with existing empirical and remote sensing data. Although nutrients are associated with sediment runoff and agriculture expansion11, we were not able to explicitly model nutrient runoff here because of the lack of spatially-explicit information on local fertilizer application rates. First, we designed land-use change scenarios that represent extreme projections of where deforestation and restoration could occur using the national land-use classification system and agriculture practices, and climate change scenarios based on historical bleaching impact in Fiji and models developed for similar latitudes in Hawaiʻi. Second, we adapted, calibrated and applied the land-sea framework, which connects the sediment and coral reef models. To measure proxies of ecological resilience, the coral reef models focused on four benthic indicators, known to respond to land-based runoff, and four fish indicators, which represent important subsistence and cultural resources13,62,65. Last, we undertook a spatial analysis to assess the impact of future scenarios on coral reefs and identify priority areas in specific watersheds, where avoiding deforestation through conservation or promoting restoration could foster coral reef resilience in the face of climate change.

Figure 2 Linked land-sea modeling framework. (a,b) Land-use and climate change scenarios were coupled with the linked land-sea framework163. (c) Land cover, topography, soil, and climate data were inputs in (d) the InVEST Sediment Delivery Ratio (SDR) model to quantify sediment export (t.yr−1). (e) The coastal discharge models used GIS distance-based dispersion models to generate the terrestrial driver grid data (t.yr−1). (f) Bathymetry and habitat maps were combined with (g) GIS-based models to derive the marine driver grid data (i.e., habitat topography, geography, exposure and complexity). (h) The coral reef predictive models were calibrated on coral reef survey data. (i) Outputs were: (1) response curves, (2) maps of benthic (% cover) and fish (kg.ha−1) indicators, and (3) a linked land-sea decision-support tool. Full size image

Human drivers scenarios

Land-use change scenarios

We considered three feasible land-use change scenarios (maintenance of current land-use, deforestation, and restoration) that represent alternative projected and extreme land-use change in Kubulau and generally in many other watersheds across tropical high islands9. The scenarios were designed based on a land-use capability (LUC) classification map developed by Fiji Department of Agriculture (Land Use Section) to inform their spatial land-use planning66. This LUC classification system is based on land development planning, soil conservation, and promotion of sustainable land-use practices. The LUC classification has eight classes: classes I–IV include land suitable for arable cultivation; classes V–VII include land suitable only for pastoral or forestry use; and class VIII is land suitable only for conservation. Each LUC class includes a range of natural limitations with respect to human uses after accounting for the physical qualities of the soil and the surrounding environment. Limitations were assessed by the Department of Agriculture based on susceptibility to erosion, steepness of slope, susceptibility to flooding, wetness, drought, salinity, depth of soil, soil texture, stoniness, structure, nutrient supply and climate, which in turn affect productivity, land-use practices, and type of land-use and associated intensity66.

In ArcGIS, we reclassified the current land-use of each pixel to future land-use based on the LUC class. Under the maintenance of the current land-use scenario, all the existing native forest is placed under conservation. The current land cover was derived from satellite imagery (see Brown et al.42 for more details). Under the deforestation scenario, we assumed that all the areas suitable for logging (LUC: V–VII) within designated current and proposed logging areas (logging concessions) were cleared and converted to plantations of pine (Pinus caribaea), an introduced timber and the principal forestry species in Fiji (described in Jupiter et al.67). Outside of logging concession areas, areas suitable for agriculture (LUC I–IV) were converted to commercial agriculture, such as monoculture of taro (Colocasia esculenta) and kava (Piper methysticum). These two root crops are the most important cash crops in this part of Vanua Levu and important drivers of deforestation across many Pacific Islands68. We assumed that existing native forest located in LUC VIII were not deforested. Under the forest restoration scenario, we assumed no loss of the current cover of native forest and all the existing pine plantations from previous logging activities were converted back to native forest, a strategy promoted by government and NGOs through payment for ecosystem services (PES) schemes to reduce soil erosion69. While most efforts include restoration with introduced species70, we assume restoration to native forest because of its conservation value71.

Climate change scenarios

Three climate change scenarios (low, moderate, and high impact in terms of levels of coral bleaching) were designed to assess the potential impact of increases in SST and atmospheric carbon dioxide (CO 2 ) on coral reefs. The coral bleaching scenarios were derived from recorded and projected coral bleaching impacts for the region72,73. An average greenhouse gas emissions scenario (A1) was assumed for the years 2000–2099 A.D. (21st century), which corresponds to a future with very rapid economic growth, global population peaks in mid-century and declines thereafter, a rapid introduction of new and more efficient technologies, and an energy system with no heavy dependence on one particular source74. Based on SST and atmospheric carbon dioxide (CO 2 ) projections, shallow-water scleractinian coral cover loss due to bleaching was estimated from a combination of growth and mortality models72. Based on a projected increase in global temperatures of 2–4 °C over the coming century with a threshold for heat stress increasing by 0.1 °C every decade, the model suggested a coral cover decline of 25% to 75% for the main Hawaiian Islands by the end of the century72,74. We compared those modeled bleaching projections to existing coral reef data collected in 2000, which monitored bleaching events impacts in Fiji73,75. The modeled predictions of coral bleaching impacts were within the same range as recorded coral bleaching impacts in Fiji73,75. Given the similar levels of impact between Hawaiʻi and Fiji, it was deemed acceptable to transfer the projected bleaching impacts developed for Hawaiʻi to Kubulau. In our low and moderate bleaching scenarios, the current coral cover for reef areas shallower than 5 m was reduced by 10% and 30%, respectively. In our high bleaching scenario, to the current coral cover between 0–5 m and 5–10 m was reduced by 30% and 10%, respectively, because deeper waters are cooler and can reduce the impact of increased SST76. Given the large uncertainties and debate surrounding coral adaptation to heat stress77, these scenarios are not to be considered quantitative forecasts of percent coral cover change for this location, but instead conservative and large-scale probability-based estimates of the relative impact of predicted increases in SST and CO 2 on coral reefs in Fiji over the next 100 years.

Linked land-sea modeling framework

The land-sea framework is made of four key components: (1) InVEST SDR, (2) coastal sediment plume models, (3) marine drivers, and (4) coral reef predictive models. Sediment export (t.yr−1) was modeled at 30 m × 30 m resolution for each watershed using the InVEST SDR, based on the land-use from each scenario, topography, soil types, and rainfall data (Fig. 2c,d). The modeled sediment export was diffused from pour points representing stream mouths (Fig. 1) into the coastal zone using GIS distance-based models to generate the terrestrial driver grid data (60 m × 60 m) (t.yr−1) (Fig. 2e). Bathymetry and habitat maps, derived from remote sensing imagery, were coupled with GIS-based models, to generate marine driver grid data at 60 m × 60 m (geography and habitat) (Fig. 2f,g). The coral reef predictive models used Boosted Regression Trees (BRT) calibrated on local coral reef survey data and generated response curves representing the relationships of each individual environmental driver to each coral reef indicator, resulting in predictive maps of the benthic (% cover) and fish (kg.ha−1) indicators (60 × 60 m) (Fig. 2h,i). Once calibrated on local data, we applied this linked land-sea framework as a decision-support tool using scenario planning to identify coral reef areas vulnerable to sediment and climate change impacts, which we traced back to identify priority watersheds for forest conservation or restoration to promote coral reef resilience (Fig. 2i).

InVEST Sediment Delivery Ratio export models

For each scenario, we applied the InVEST SDR model59 to quantify the sediment export (t.yr−1) to the coast by watershed due to soil loss on hillslopes from overland erosion59,60 (Fig. 1c). SDR is spatially-explicit and operates at the resolution of the digital elevation model (DEM) input (30 m). For each cell, the model computes the amount of eroded sediment using the revised universal soil loss equation (RUSLE), and a sediment delivery ratio (SDR) to estimate the proportion of soil eroded on a given area that will travel to the stream mouth at the shoreline (see Hamel et al.59 for full details on the model). This approach was initially proposed by Borselli et al.78, which relies on modeling sediment transport throughout the landscape based on local topography, and therefore does not require hydrological modeling to determine the sediment ratio exported to the shoreline.

First, we estimated the overland gross erosion per cell using the empirical Revised Universal Soil Loss Equation (RUSLE) (see Equation 1)79:

$${Soil}\,{loss}={R}\times {K}\times {LS}\times {C}\times P\,$$ (1)

where R = rainfall erosivity (MJ.mm.ha−1.hr −1), K = the rate of soil loss per rainfall erosion index unit, known as soil erodibility (ton-ha-hrs.MJ−1.ha−1.mm−1), LS = slope-length and gradient factor (derived from the DEM [30 × 30 m]59), C = a vegetation cover (C-factor) and P = management practice effectiveness (P-factor) (see Table S3 for more details on parameters). The map of rainfall erosivity was derived from monthly rainfall averages80 at a 100 × 100 m (P) and converted to erosivity using the Bols method applied in Indonesia (see Equation 2)81:

$${R}=\frac{2.5\times {{P}}^{2}}{100(0.073{P}+0.73)}$$ (2)

Soil erodibility (K), the rate of soil loss per rainfall erosion index unit82,83, was derived from the New Zealand Soil Survey dataset84, and used a value of K of 0.002 ton.ha.hr.MJ−1.ha−1.mm−1 to fill in missing values that were not available in the tables. C-factors derived from existing literature and studies conducted in similar regions were assigned to each land cover/use type (see Table S3 for more details on parameters)82,85,86,87,88. Management practice effectiveness (P factor) was not considered in this model due to lack of data59.

Then, the SDR model estimated the proportion of soil eroded on a given area that travelled to the stream mouth at the shoreline. The SDR model is based on the concept of hydrological connectivity to estimate sediment retention and export to the shoreline (see Borselli et al.78 for more details). First, the SDR computes a connectivity index (IC i ) for each pixel i based on the upslope area and downslope flow path78. A streamflow accumulation threshold was set to define streams based on the DEM59. Given the lack of empirical data for the region, the connectivity of the model was verified by comparing predicted stream outputs to available stream maps. Sub-watersheds were created using the Basins function in ArcGIS and pour points at the shoreline were edited for accuracy in comparison to satellite imagery and bathymetry data89,90. Then, a sediment delivery ratio was derived for each pixel i based on the connectivity index59. The SDR model parameters include an IC 0 , a Borselli k-factor, and a maximum allowable SDR that define the shape of the relationship between the SDR and the connectivity index (IC i )59. The calibration parameters IC 0 and the k-factor were set to 0.5 and 2.0, respectively, and the maximum allowable SDR was set to 0.859 (see Sharp et al.58 for more details on effects of parameterization). This approach was selected since it requires a minimal number of parameters and is spatially-explicit. We note that the models have yet to be quantitatively validated against local datasets and were parameterized with values from other regions, which can differ in terms of climate and soil conditions91.

Coastal sediment plume models

To link the sediment and coral reef models, we modeled the impact of sediment runoff on coastal water quality. In order to do so, we generated a water quality map (60 × 60 m) representing the total sediment load (represented by TSS) from all the watersheds discharging into Kubulau coastal waters, for each land-use scenario. First, the modeled sediment export from each watershed was diffused into coastal waters by adapting a previously developed dispersal plume model in ArcGIS to represent the point source nature of stream discharge in local coastal water conditions53,57,92 (Fig. 2d). To do so, we created a cost-path surface (c) that quantifies the least accumulative cost-distance (impedance) of moving planimetrically through each cell from each pour point using a composite of three marine drivers known to affect diffusion: depth (m), distance to stream mouth (m), and wind exposure (degree) (see ‘Marine drivers models’ section below for more details)11,57,93. Then, the spread of sediment into coastal waters from each pour point was modeled using a decay function, which assigned a portion of the remaining quantity from the previous cell in all adjacent cells, based on the cost-path surface until a maximum distance of 3 km from the shoreline was reached53,57 (see Equation 3). This threshold was based on the distance between river mouths measured in ArcGIS and locations where sediment impacts on coral reefs have been recorded in the past94.

$${{S}}_{{i}}={{s}}_{{p}}\times {{e}}^{-{{c}}^{2}/{{D}}_{{c}}}$$ (3)

where S i = cell value for dispersed sediment (t.yr−1) for watershed i, S p = Sediment load (t.yr−1) at each watershed pour point (obtained from summarizing the total sediment export per watershed), c = cost-path surface (unitless), D c = cost-path surface threshold distance from the shore for each decayed sediment plume per watershed (equivalent to 3 km from the shoreline). Last, we summed all the individual watershed sediment plume gridded maps in ArcGIS to obtain the total sediment load (represented by TSS) per land-use scenario for each pixel of coral reef area. This approach to modeling coastal sediment discharge is diffusive, and thus allows for wrap around coastal features, but does not account for nearshore advection that acts to push suspended sediment in specific directions53. We used these diffusive models to derive conservative estimates of sediment plumes, since the nearshore circulation patterns were unknown for our study site.

Coral reef indicators

To measure coral reef resilience, we considered the percent cover of four benthic groups (crustose coralline algae [CCA], scleractinian corals, turf algae and macroalgae) and the biomass of four fish groups (kg.ha−1) based on their ecological roles and cultural importance to local communities62 (see Table S1 for more details). The benthic indicators are known to respond to changes in land-based runoff, which influence the distribution of fish taxa11,40, and therefore support aspects of coral reef ecological resilience13,65. Fish taxa identified as important for subsistence and cultural practices by the local communities were modeled according to their ecological role: (1) browsers, (2) grazers, (3) scrapers, and (4) predators (see Tables S1 and S2 for more details)62,65. We derived percent cover of the benthic indicators (%) and biomass of the fish indicators (kg.ha−1) from reef survey data collected by the Wildlife Conservation Society (WCS). The field dataset comprised 163 survey locations randomly stratified by depth (deep [12–15 m], shallow [5–8 m], top [0.5–2 m]), habitat (forereef and backreef areas), and management (open or closed to fishing) (Fig. 1), collected over three sampling periods (2009–2010) (see64,95 for more details).

Marine driver models

The marine driver grid maps (60 × 60 m) were derived from bathymetry96 and benthic habitat95 maps for the site using GIS-based tools (Fig. 2f,g)97,98. Based on existing literature and local community input, local geography (depth and distance to shore) and habitat characteristics (topography, complexity, and exposure) were identified as important drivers of the modeled coral reef benthic and fish indicators (see Table S4 or see Stamoulis & Delevaux et al.97 for more details on processing methods). Depth and distance from shore were used to account for variation arising from spatial location. Depth was derived from bathymetry at 4 m resolution96 using passive remote sensing techniques, and distance from shore was calculated as the Euclidean distance (m) from the coastline, using the national coastline map in ArcGIS99. Three types of habitat drivers, representing direct and indirect effects of seafloor topography on benthic and fish communities were also derived from this bathymetry data96: (1) habitat topography, (2) habitat complexity, and (3) habitat exposure. Habitat topography metrics, represented by Bathymetric Position Index (BPI) and slope, described the position of the reef relative to the surrounding area. These metrics were computed using the Benthic Terrain Modeler tool for two neighborhood sizes (60 and 240 m radii) to capture habitat topography at different spatial scales100,101,102. Habitat complexity metrics that describe fine-scale topographic structure, represented by slope of slope, planar curvature, and profile curvature were computed using the DEM Surface and Curvature Tools in ArcGIS96,98. Habitat exposure metrics were used to characterize the direct and indirect effects of water flow due to fine-scale seafloor topography and directionality. These metrics were derived by computing seafloor aspect, the steepest downslope direction of the seafloor measured in degrees, using the Aspect tool in ArcGIS96,98. We estimated exposure to winds by calculating the circular mean and standard deviation of aspect and converted the circular mean into measures of northness and eastness using the sine and cosine functions, respectively, from the Spatial Analyst toolbox in ArcGIS96,98. Four types of habitat connectivity metrics, representing direct and indirect effects of habitat composition and fragmentation on benthic and fish communities, were derived from the benthic habitat map, with a 10 m minimum mapping unit in FRAGSTATS software: (1) contiguity, (2) fractal dimension, (3) proximity, and (4) Shannon diversity index95,103. We note that it is not necessary to generate all these drivers in order to apply this method elsewhere, it depends on whether a bathymetry and/or habitat map is available.

Spatial predictive coral reef models

We used BRT to build the coral reef models (Fig. 2h)104. Tree-based models are effective at modeling nonlinearities, discontinuities (threshold effects), and interactions between variables, which is well suited for the analysis of complex ecological data105,106,107. BRT models can accommodate many types of response variables, including presence/absence, count, diversity, and abundance data108. Since the coral reef indicators were all continuous variables, the response variables were modeled using a Gaussian (normal) distribution, and appropriate data transformations (square root for benthic indicators and fourth root for fish biomass) were applied to improve the normality of the response variable distributions. The highly correlated (r > 0.7) environmental drivers were removed from the coral reef models97. We calibrated the BRT models on coral reef data to determine the most influential drivers (among the simultaneously tested predictors) and estimate the underlying relationships between the modeled indicators and the key drivers using response curves108,109. The values of the terrestrial and marine drivers’ grid maps were sampled using bilinear interpolation at the location of each reef survey (start of the transect) in ArcGIS. This approach takes a weighted average of the 4 nearest cell values, thereby accounting for the relative position of the reef surveys on the predictor grids and their different native spatial scales.

A BRT model was independently developed for each coral reef indicator to determine the relative influence of terrestrial and marine drivers using the values of the coral reef indicators and interpolated terrestrial and marine driver values at the reef survey locations. First, we calibrated the benthic indicators’ models as a function of the terrestrial and marine drivers. Then, we calibrated the fish indicators’ models as a function of the terrestrial and marine drivers, as well as the empirical abundance of benthic groups as additional predictors in the models for the fish groups. The calibration process used an internal ten-fold cross-validation to maximize the model fit and determine the optimal combinations of four parameters: (1) learning rate (lr), (2) tree complexity (tc), (3) bag fraction (bag), and (4) the maximum number of trees (see Elith et al.108 for more details). We used the percent deviance explained (PDE) and internal ten-fold cross validation PDE (CV PDE) as performance measures of the model optimum. The optimal models explained the most variation in the response variables (i.e., greatest CV PDE) and were selected as the best and final models. The model calibration was conducted in R software using the gbm package108,110,111. Spatial autocorrelation of the response variables was tested using Moran’s I Index for both the raw values and the ecological model residuals112.

Using the calibrated coral reef models, we predicted and mapped the distribution of each coral reef indicator on a cell-by-cell basis using the values of the terrestrial and marine drivers at each grid cell across the coral reef model domain. Predictive maps were generated for each indicator under present conditions and future scenarios (land-use change, climate change, and combined scenarios), described in more detail below (Fig. 2a,b)72,73. The coral reef predictive maps were generated at 60 × 60 m to account for the dimensions of the reef survey methods (i.e., 25 m transects), the transect bearings, and the positional accuracy of global positioning system used to navigate to them in the field113,114. The boundaries of the coral reef model domains comprised the qoliqoli boundaries to capture the spatial extent of this management unit and the offshore boundary corresponded to the maximum surveyed depth (i.e., 22 m) (Fig. 1b). First, we spatially predicted each benthic indicator as a function of its key drivers. Then, we spatially predicted the fish indicators as a function of their key drivers, including the predicted distribution of the benthic indicators.

In order to evaluate the quality of the coral reef model predictions, we compared the measured and predicted values of the coral reef indicators under present conditions. The values of the interpolated predictions and surveyed coral reef indicators at these locations were compared with a linear regression (R2 and p-value). The predicted values of the benthic and fish grid maps were sampled using bilinear interpolation at the location of each reef survey (start of the transect) in ArcGIS, thereby accounting for the relative position of the reef surveys on the predicted grids. Then, we overlaid the predicted maps with the survey point values for each indicator using the same color ramp scale for the legend to enable visual comparison. The spatial predictions were performed in the R software using the dismo and raster packages111,115,116.

Coral reef scenario impact assessment

We identified coral reef areas that could be impacted by land-use and/or climate change scenarios. To do so, we calculated differences between predictions of coral reef indicators under the land-use change, climate change, and combined scenarios, compared to present conditions. We computed the significance of the pairwise differences per grid cell for each coral reef indicator relative to the mean and variance of all differences across the coral reef model domain using the SigDiff function from the R package SDMTools117. The grid cells representing significant differences (α = 0.10) were reclassified to indicate where predictions were significantly different than present conditions under each scenario34. For the combined scenarios, we overlaid the coral reef areas of significant differences under land-use and climate change scenarios and delineated areas of overlap, where both drivers could potentially interact. These areas were combined into a single map to display the spatial pattern of potential impact per scenario. Finally, the areas of significant differences for each coral reef indicator were used to quantify the relative changes in benthic habitat and fish biomass within those areas. For areas where coral bleaching and sediment runoff overlapped, we summed the changes in abundance of each coral reef indicator. Although interactions between sediment runoff and bleaching were not explicitly modeled, this approach allowed us to delineate coral reef areas where cumulative impacts and interactions could occur. High turbidity and suspended sediment can mitigate elevated SST and associated bleaching impacts (i.e., antagonistic effect) but these interactions remain poorly understood16,17. Due to lack of data, nutrient runoff was not explicitly modeled in this study, but nutrients are known to bind and travel with sediment, thereby potentially contributing to lack of recovery from bleaching through promoting algae growth (i.e., synergistic effect)14,15. However, given the limited understanding of these complex processes and the lack of data to represent nutrient runoff, we assumed that these effects were additive. Therefore, we applied the precautionary principle and considered these coral reef areas as vulnerable to cumulative impacts118.

Spatial prioritization of conservation

In order to locate the most effective areas to prioritize forest conservation (i.e., prevent deforestation) and restoration to foster coral reef resilience, we linked the coral reef areas vulnerable to sediment runoff from land areas within upstream watersheds that contributed the major portion of the sediment load to those areas. Using the individual watershed sediment plume maps (60 m × 60 m) from the coastal plume models, we identified the watersheds contributing the majority of sediment runoff to coral reef areas likely to change under land-use and/or climate change scenarios, compared to present conditions. To do so, we linked the coral reef areas showing significant difference to the watersheds upstream, which contributed the majority (>90%) of the total sediment load (represented by TSS) at those areas in R and ArcGIS. The linked coral reef areas and watersheds were combined into a single map to display the land-sea linkages for the land-use and climate change scenarios. Then, using the sediment export maps (30 m x 30 m) from the SDR models, we identified land areas contributing the most sediment runoff to coral reef areas likely to change under land-use and/or climate change scenarios, compared to present conditions. For climate change scenarios, we identified priority land areas in the linked watersheds contributing a large portion (>66%) of sediment export to downstream coral reefs under current land-use. For the land-use change scenarios, we calculated the significant differences in sediment export in the linked watersheds per grid cell, compared to present conditions using the SigDiff function117. For each scenario, the grid cells representing significant differences were reclassified to indicate where sediment export was significantly different from present conditions (α = 0.10). In ArcGIS, we created a 100 m buffer around those priority land areas based on conservation and logging best management practices60,119. To visually represent those results, the linked watersheds were combined into a single map to display where conservation and/or restoration actions could foster coral reef resilience based on the land-use and climate change scenarios and overlaid the priority land areas for each scenario.

Data availability statement

All the data produced in this study is available from the data repository fisgshare: https://doi.org/10.6084/m9.figshare.6823508.