Abiotic and Human Drivers of Reef Habitat Complexity Throughout the Main Hawaiian Islands

Reef rugosity, a metric of three-dimensional habitat complexity, is a central determinant of reef condition and multi-trophic occupancy including corals, fishes and invertebrates. As a result, spatially explicit information on reef rugosity is needed for conservation and management activities ranging from fisheries to coral protection and restoration. Across archipelagos comprising islands of varying geologic stage and age, rugosity naturally varies at multiple spatial scales based on island emergence, subsidence, and erosion. Reef rugosity may also be changing due to human impacts on corals such as marine heatwaves and nearshore coastal development. Using a new high-resolution, large-area mapping technique based on airborne imaging spectroscopy, we mapped the rugosity of reefs to 22 m depth throughout the eight Main Hawaiian Islands. We quantified inter- and intra-island variation in reef rugosity at fine (2 m) and coarse (6 m) spatial resolutions, and tested potential abiotic and human drivers of rugosity patterns. We found that water depth and reef slope remain the dominant drivers of fine- and coarse-scale rugosity, but nearshore development is a secondary driver of rugosity. Our results and maps can be used by fisheries management and reef conservations to track geologic versus human impacts on reefs, design effective marine managed areas, and execute activities to improve reef resilience.


INTRODUCTION
Habitat complexity is the three-dimensional (3D) structure of the physical environment with which organisms interact and is a key influencer of the distribution of marine biota. Habitat complexity can affect a variety of ecological functions such as by providing prey protection from predators, surfaces for coral settlement and growth, and different microhabitats to support a variety of species. Habitat complexity interacts with other environmental factors such as light, temperature and wave action to generate a wide range of ecological conditions (Sebens, 1991).
On coral reefs, 3D habitat complexity is a determinant of reef fish assemblages, where fish abundance, biomass, and richness are often positively correlated with complexity (Cinner et al., 2009;Graham et al., 2009). Yet the direction and strength of this relationship does vary, and studies can also show mixed connectivity between reef complexity and fish stocks (Jennings et al., 1996;Harborne et al., 2012). Coral cover is also often positively correlated with reef complexity (Alvarez-Filip et al., 2009;Graham and Nash, 2012). Correspondingly, negative relationships between complexity and algal cover often exist, perhaps indirectly reflecting the role of reef complexity in predicting reef fish and resulting herbivory (Graham and Nash, 2012). However, 3D habitat complexity is not only driven by coral cover, but also by benthic geomorphology at a range of scales from small boulder and erosion-deposition zones to large subsurface geologic structures.
The Main Hawaiian Islands (MHIs) are an important casein-point, varying widely in geologic age from less than a few years old on Hawai'i Island to more than six million years old on Ni'ihau (Neall and Trewick, 2008). Island age is accompanied by stage of accretion and subsidence, processes that generate enormous inter-and intra-island variation in reef extent associated with benthic substrate availability (Asner et al., 2020b). Younger islands such as Hawai'i contain vast fringing reefs dominated by rock substrates (i.e., lava beds, massive boulders, rock fields) that define much of the reef rugosity, but with embayments and coves co-dominated by geology and coral growth. In contrast, older islands such as Kaua'i and O'ahu contain reefs defined by geologic subsidence and erosional surfaces, with broad expanses of both sandy bottom and calcareous reef substrate.
Overlain on this geologic template, water motion and light penetration decrease with increasing depth, and wave-sheltered areas show the greatest Holocene reef accretion (Grigg, 1998). Generally, optimal reef growth occurs between 10 and 20 m depth in the MHI, reflecting a trade-off between wave-induced stress and decreased light availability (Storlazzi et al., 2005). In areas of greatest wave power, turf algae and crustose coralline algae dominate due to more wave-tolerant, flatter morphologies. On the other hand, low wave power regions are dominated by hard corals and macroalgae with structures more vulnerable to dislodgement (Engelen et al., 2005;Madin et al., 2014;Gove et al., 2015).
Reef complexity varies over multiple spatial scales (Mumby et al., 2004). Measurements of complexity aim to quantify the vertical variation of the benthic surface in relation to the two-dimensional (2D) surface of the same area, and therefore complexity metrics are dependent upon the resolution of the mapping technique. Larger grid or pixel sizes of bathymetric maps result in complexity values describing broader geological features in comparison to smaller pixel sizes that represent specific coral colonies, rock outcrops and other fine-scale features. Presently, the spatial scale at which biophysical drivers influence fish or invertebrate assemblages is poorly known (Purkis et al., 2008;Aston et al., 2019), and a multi-scale assessment of reef complexity would improve our understanding of habitat across a range of scales relevant to a wide range of ecological processes (Harborne et al., 2006;Williams et al., 2015).
Over the past several decades, there has been a rapid increase in the diversity of methods used to quantify 3D habitat complexity, where various techniques are aimed at quantifying complexity at specific scales and spatial resolutions. Studies focusing on colony-level measures most commonly use a chain-and-tape method to estimate rugosity, an index of habitat complexity. Rugosity is estimated by calculating the ratio of chain length laid across the bottom reef profile to the linear distance of the transect (Risk, 1972). Complexity has also been estimated visually by scoring a range of structural variables (Gratwicke and Speight, 2005). Modern methods of assessing habitat complexity include 3D computer models, which are generated from images via structure-from-motion photogrammetry, allowing quantification of complexity metrics such as surface rugosity among fractal dimension, total surface area and surface height (Burns et al., 2015;Figueira et al., 2015). Field techniques, however, can be time-consuming and typically operate over limited areas (<0.25 ha; Lechene et al., 2019).
Remote sensing methods-primarily acoustic imaging and light detection and ranging (LiDAR)-have allowed estimations of reef complexity over areas of much greater extent. More than 50 years ago, sound navigation and ranging (SONAR) revolutionized reef mapping over large areas by providing moderate resolution information (Moravec and Elfes, 1985). Subsequent LiDAR-derived rugosity measures of reef at 4-5 m spatial resolution have been correlated with in situ measures of reef rugosity (Wedding et al., 2008). However, measurements of habitat complexity coarser than 10 m did not show significant relationships with in situ rugosity, but they still provided relevant information at broader geographic scales (Wedding et al., 2008). Similarly, reef rugosity estimated though acoustic methods showed relationships with fish at kernel sizes less than 40 m, with the most significant relationships with species richness at a resolution finer than 8 m (Purkis et al., 2008). More recently, airborne imaging spectroscopy has opened up access to large-area bathymetric and reef rugosity mapping at resolutions of 0.4 m (Asner et al., 2020a).
Throughout the Hawaiian Islands, it is not known how environmental and physical factors compare with human influences in shaping 3D habitat complexity (Wedding et al., 2018). Understanding controls over spatial patterns of complexity is relevant to creating effective marine managed areas, where planning can benefit fish assemblages, protect essential fish habit and enhance recruitment (Grober-Dunsmore et al., 2007). This is especially important for coastal communities that heavily rely on ecosystem services from reef fish populations. Furthermore, protecting locations with greater reef complexity should be a management priority as they have increased resilience to coral bleaching and enhanced recovery through larger fish populations (Graham and Nash, 2012;Rogers et al., 2014). Considering the wealth of research linking reef complexity and reef fish, rapid methods that can document rugosity at broader geographic scales are of great value. Here we use a new airborne mapping technique to report on how 3D habitat complexity varies at fine and coarse spatial resolutions within and between each Main Hawaiian Island. We undertook this study without a specific focus on live corals, which was largely covered in Asner et al. (2020b). In that earlier paper, we determined the geographic distribution of live corals, but we were unable to associate live coral cover broadly with rugosity due to extremes in island age, geologic stage, and erosional processes. As a result, in the present submission, we determined the geographic distribution of fineand coarse-scale rugosity surfaces in order to assess on the relative importance of multiple abiotic environmental as well as human factors potentially affecting reef complexity at intra-and inter-island scales.

Bathymetric Mapping
Using data collected by the Global Airborne Observatory (Asner et al., 2012), we mapped benthic depth at 2 m spatial resolution throughout the eight MHI. We selected 2 m mapping resolution as a practical trade-off between time and cost of mapping the vast area of the MHI, while simultaneously achieving sufficient resolution to resolve fine-scale variation in reef structure that incorporates larger coral colonies and basaltic rock outcrops (see Asner et al., 2020a). The GAO collected data from multiple coaligned instruments, two of which were used for bathymetric and rugosity mapping: a high-fidelity visible-toshortwave infrared (VSWIR) imaging spectrometer and a dualbeam LiDAR scanner (Asner et al., 2012). Data from the spectrometer were used to model benthic depth using the deep learning methodology described, tested and validated in Asner et al. (2020a).
The VSWIR spectrometer and LiDAR data were collected between January 2 and February 4, 2019. To maximize data consistency, daily airborne operations were performed from 0830 to 1100 local time. Collection location could change during each flight day and was actively managed based on need, cloud cover, and windspeeds to provide both the most efficient use of time and the best conditions for spectroscopic seafloor measurements. During flight, instrument settings were set for the planned nominal flight altitude of 2 km above the sea surface. Flightlines were spaced to achieve 50% overlap in VSWIR spectrometer coverage. Aircraft ground speed was 130-140 kt. LiDAR pulse frequency was set to 200 kHz (100 kHz per laser) and scan frequency was 34 Hz with a field of view of 38 • , allowing 2 • of buffer on each side of the spectrometer field of view of 34 • , achieving a nominal pulse density of more than 4 pulses m −2 . The radiance data from the spectrometer are collected in 427 spectral channels covering the 350-2500 nm wavelength range in 5 nm increments. Using a modified version of the ATREM model (Gao and Goetz, 1990;Thompson et al., 2017), we retrieved ocean surface reflectance from the at-sensor radiance data. The reflectance data for each flight line were passed through the model described in Asner et al. (2020a) to retrieve estimated depth in meters after masking out problematic regions (clouds, waves, and regions of excessive solar glint) from the data. Because water absorbs nearly all infrared light, land surface objects were filtered from analysis using a reflectance threshold at the 890-910 nm wavelength.
Orthorectification of the spectrometer data was a multi-step process. During flight GPS timing data were collected during flight by both instruments, and this timing location was used to link the data to a precise flight trajectory built from corrected and interpolated data from the onboard positioning and orientation system (POS). With a LiDAR-derived digital surface model (DSM) and a known position and orientation of the sensor at each sample, the 3-dimensional position of each spectrometer pixel was ray-traced to the sea surface level. However, to accurately map the location of each pixel on the ocean floor, we needed to account for estimated depth (d) and the refraction of light at the ocean surface. For each pixel, observation angle from zenith at the water surface, φ a , is known and with standard refractive indices of 1.33000 for water and 1.00029 for air, Snell's refraction formula can be used to compute the angle after crossing the water surface boundary as φ b = sin −1 (0.75210 sin φ a ). This angle can be traced down to the ocean floor to retrieve the 3-dimensional ocean floor location for each pixel.
We took advantage of the 50% overlap in flight lines to build smooth depth mosaics by reducing sharp transitions between flight lines. For each group of flight lines making up a single region of coastline, the process for mosaicking followed these steps: (1) Record pixel id, flight line id, depth, and observation zenith angle for up to 5000 (randomly sampled, if more than 5000 are available) pixels that have overlapping flightlines. This was done for each coastal region for each of Ni'ihau (one region), Kaua'i (four regions), O'ahu (four regions), Moloka'i, Lana'i, Kaho'olawe (each one region), Maui (four regions), and Hawai'i (11 regions). These individual regions were then hand-mosaicked into whole-island maps of depth, with the exception of Hawai'i which was broken into four quadrants due to size limitations.

Multi-Resolution Rugosity
We computed island-wide maps of rugosity at two resolutions using the "planar" rugosity metric (Jenness, 2004) on the GAO benthic depth maps (Figure 1). Missing data of less than two pixels in width were filled using an inverse-distance weighted average of the three nearest neighboring pixels. Finescale rugosity was computed using a 3 × 3 pixel (6 × 6 m) moving window on the original 2 m resolution bathymetric maps (Figure 1b). This resolution seeks high frequency changes on the seafloor arising from coral colonies, rocks, and other bottom features that generate local habitat variability. Coarse-scale rugosity was computed by first down-sampling the 2 m depth maps to 6 m resolution using a mean filter. The planar rugosity metric was then computed using a 9 × 9 pixel (54.0 × 54.0 m) moving window on the 6 m depth maps. Coarse-scale rugosity is responsive to variations in larger terrain features resulting from geologic, reef-scale accretion and subsidence processes. Previous testing at other resolutions indicated that a 9 × 9 pixel coarsescale rugosity product gives the best detail for geologic features while removing fine-scale "noise" from corals and rock boulders (Figure 1c). These two resolutions of rugosity were computed for each of the individual island depth maps. Because the distribution of rugosity (r) is heavily skewed to the right because of occasional noise pixels, and the units of the output are not meaningful, we transformed raw rugosity values to have an approximate uniform [0,1] distribution T by sorting the n individual values across the map and assigning the value of transformed rugosity Thus, for the island-scale fine and coarse rugosity maps, low transformed rugosity values will be near-zero, mid-range values will be near 0.5, and extreme values will never be more than 1.0.

Geospatial Analyses
We assessed the spatial patterning of fine and coarse rugosity both between and within the MHIs using empirical variograms. To overcome computer memory and time limits, we broke the rugosity maps for each island into a grid of 1 × 1 km cells, made up of 500 × 500 pixels for fine rugosity and 167 × 167 pixels for coarse rugosity. Because of the additive nature of the variogram variances, we could compute one for each island using samplesize weighted averages of the individual lag variances across grid cells into a single variogram. This methodology provided the additional benefit of generating information about the variation in spatial patterning within each island.
For each grid cell, the associated square region was extracted from the rugosity map along with a buffer of 500 m (250 pixels for fine-scale maps and 83 pixels for coarse-scale maps) on all sides. The empirical variogram for each grid cell was computed with the GeoStats package (Hoffimann, 2018) in the Julia programming language (Bezanson et al., 2012), using 40 lag steps and a maximum lag of 500 m. A variogram model was fit to the empirical variogram for each grid cell, with the fitting software iterating over eight model forms to find the optimal form as well as the associated optimal estimates for the model parameters: sill, nugget, and range parameters. By storing the sill, nugget, and range value for each parcel, we were able to aggregate these into a 1 × 1 km resolution map of each parameter for each island. The island-wide empirical variogram was computed as a weighted average of the variograms from all grid cells, where the weight for each cell was the number of valid rugosity observations separated by the given lag distance within that cell.

Land-Sea Driver Modeling
We used multiple land-sea environmental maps available for the MHI to assess the relative importance of environmental and human factors that may contribute to the mapped distribution of fine-and coarse-scale rugosity ( Table 1). For example, reef depth, slope, and distance to coastline vary with island geologic stage, thereby affecting erosion, subsidence, and other factors known to shape larger benthic structures (Fletcher et al., 2008;Gove et al., 2015). Abiotic drivers also play a role in shaping rugosity, either via physical impacts of waves and wind action (Li et al., 2016), or via the impacts of temperature on calcareous organisms like corals that contribute to rugosity (Ignatov, 2010;Wedding et al., 2018). Finally, nearshore development may have an effect on reef rugosity via sedimentation and/or removal of reef structures (Center for International Earth Science Information Network-Ciesin, 2018). We checked for driver variable covariance and utilized a combination of drivers each with less than 50% correlation.
While the GAO fine-and coarse-scale rugosity maps have a resolution of 2 and 6 m, respectively, many of the environmental maps in Table 1 are provided only at resolutions of 30 m or coarser. The GAO rugosity maps were therefore coarsened to 30 m resolution using a mean filter, because we determined this resolution to sufficiently balance the finer and coarser datasets in this analysis. Similarly, all input factors were resampled to match this resolution, using cubic spline interpolation for maps requiring upscaling and a mean filter for maps requiring downscaling. These input factor maps included the GAO bathymetry map used to compute rugosity, as well as slope computed from the 30 m depth map. We sought to analyze the relative importance and influence of each of these factors on each of the two scales of rugosity at two levels: once separately for each individual island and once for all islands combined. By separating these analyses, we assessed whether drivers of rugosity differ by island while simultaneously revealing general patterns that apply to all islands. We randomly selected 80,000 pixels from the 30 m resampled maps for each individual island and, separately, randomly selected 100,000 pixels across all islands. Modeling was carried out using a Random Forest Machine Learning (RFML) approach (Breiman, 2001) with the Scikit-Learn python package version 0.22.1 (Pedregosa et al., 2011). Optimal settings for the RFML metaparameters defining the number of trees in the forest ("nest") were discovered through a grid search approach separately for each island ( Table 2). We manually specified several metaparameters: four as the maximum tree depth, 1 as the minimum number of samples per leaf, two as the minimum number of samples per split, and four as the maximum number of features scanned at each split. All other meta parameters were left at the default settings. The final model was a descriptive machine learning model, fit to better understand the system rather than be applied to any new data. Because the risk of overfitting is negligible in these circumstances, the idea of parsimony is not applicable and no model simplification procedure was required. With the RF model-fitting procedure, unimportant factors are simply not selected as splitting criteria for many of the node splits that make up each component regression tree. Thus, these factors will have little influence over model capacity to fit the data, and The R-squared value after fitting the model to the full dataset is shown as baseline R 2 used in the permutation importance procedure.
this lack of influence will show as low values in any following assessments of individual factor importance.
We assessed the importance of each variable using a permutation-based, correlation coefficient (R 2 )-reduction metric. To reduce the complexity of the permutation procedure, the baseline R 2 value used for this analysis was first obtained by refitting the final model for each island and rugosity scale using the full dataset (no train/test split) and comparing the predictions of this full model against the map value of rugosity for all samples in the dataset (Table 2). Next, the following permutation importance analysis was performed for each of the island-level models and the combined models at each scale of rugosity: For each input variable and for each of five iterations, the values for this variable were randomly shuffled, keeping values of all other variables intact. During each iteration, predictions were again obtained using the model on the full dataset containing the permuted variable, and we retained the difference between the original R 2 and the R 2 computed from this permutation. The five difference values were averaged for each variable to get a single importance value, where larger positive values indicate greater reduction in R 2 and, equivalently, greater variable importance.
In addition to the importance of each variable, we measured the marginal trends between each of the input variables and modeled rugosity values. These trends help understand why and how the individual input variables are important, a difficult task considering the high dimensionality of these systems. As with the importance analysis, this analysis was done for each input variable for each of the models in Table 2. For each model trained with the full dataset, the predicted rugosity values for the entire dataset were computed and stored. Next, the full range of each input variable was split into 100 equally sized bins based on percentiles of the input variable. For each of these partitions, the mean model-predicted rugosity value was computed along with the median value of the input variable within the partition. Plotting the predicted rugosity against median input variable values provides a view of rugosity response relative to the variable, inclusive of all correlations with other variables in the model. Where importance can inform as to how strongly rugosity shifts with a variable, the marginal plots can show whether the relationship is positive or negative and linear or curvilinear.

Inter-Island Reef Rugosity
Across the MHI, the total shallow reef area mapped was 98,344 ha, covering about 95% of all known reefs to a depth of 22 m Figure 2 and Appendix Figure A1). On a whole-island basis, Hawai'i and Kaho'olawe contained reefs with the greatest 3D habitat complexity, with relative rugosity values in the 0.8-0.9 range ( Table 3). O'ahu had the lowest reef complexity overall, at 0.55 or nearly half that of the structurally most complex islandscale reefs. Fine-and coarse-scale rugosity roughly tracked one another at the island level.

Geospatial Variation in Reef Complexity
Despite strong inter-island differences in 3D habitat complexity, any given local reef area displayed a wide range of rugosity values (Figure 3 and Appendix Figure A2). Sandy areas with low rugosity are shown in blue, while highly complex coral and rock reef patches are expressed in yellow to red colors. Relative to shoreline, the location and spatial distribution of high-rugosity benthic surfaces was highly variable, although rugosity can be seen to broadly increase with distance to shore (Figure 3).
At the individual island level, both fine-and coarse-scale rugosity values were non-normally distributed (Figure 4). Hawai'i and Kaho'olawe were highly skewed in the positive direction, indicating the widespread presence of structurally complex benthic surfaces. Field observations indicated that these FIGURE 3 | Zoom images from maps shown in Figure 2 demonstrating high spatial frequency variation and pattern in the fine-scale rugosity. These seven example reefs encompass the range of conditions among high-rugosity sites. Coarse-scale rugosity is provided in Appendix Figure A2. areas are dominated by rock and/or coral. By contrast, O'ahu was largely comprised of low-rugosity substrate, with few rocks and coral colonies, and large swaths of sand-and algal-covered surfaces (see Asner et al., 2020b). These contrasting distributions are both an expression of island age, size, and stage as well as how much live coral can be found in each reef ecosystem. Notably, Moloka'i displayed a bi-modal rugosity distribution, with high values dominating the windward north coast and lower values found along the leeward south coast (Figure 2). Reef rugosity showed a highly variable spatial arrangement at intra-and inter-island scales (Figure 5 and Appendix Figure 3). Variogram-range maps reveal that, on geologically older islands such as Ni'ihau, Kaua'i and O'ahu, reef rugosity mostly varied at low spatial frequency, often more than 100 m, indicated by high variogram-range values (Figure 5a and Appendix Figure 3a). However, this background condition was punctuated by reef locations with relatively high spatial frequency changes in rugosity. By contrast, on geologically younger islands such as Hawai'i, Lana'i and Maui, low variogram-range values indicated high-frequency variability in rugosity over large tracks of reef. Variogram-sills represent the baseline variation in rugosity between areas on an island that are far enough apart to lack any shared local environmental controls. Again, geologically older islands were largely uniform (Figure 5b and Appendix Figure 3b), whereas younger islands were far more variable. Overall, based on island-level semi-variograms, Lana'i and Hawai'i were among the most variable in terms of fine-or coarse-scale rugosity (Figure 5c and Appendix Figure 3c).

Drivers of Habitat Complexity
Machine learning analyses revealed a heirarchical set of contributors to the spatial variability of fine-and coarse-scale rugosity throughout the MHI (Figure 6). Fine and coarse models accounted for 86 and 73%, respectively, of the overall mapped geospatial variation, indicating that the driver variables used in the model were appropriate in assessing rugosity at both resolutions. Water depth and reef slope exerted the strongest influences on the spatial variation in rugosity as evidenced in model permutation reduction R 2 -values of up to 0.88 and 0.69, respectively, among all islands combined for fine-and/or coarsescale rugosity (Figure 6). Secondarily, nearshore development and distance from coast exerted a detectable effect (R 2 up to 0.31). Other variables such as sea surface temperature, wind speed, total land effluent, and substrate age accounted for far less variation in rugosity (R 2 < 0.11).
There was a wide inter-island range of influence of each major driver on 3D habitat complexity. First, depth was far more important than reef slope in driving fine-scale rugosity, but depth alone ranged in relative importance from low on Hawai'i to extremely high on O'ahu (Figure 6a). In contrast, reef slope effect on fine-scale rugosity was minimal on Ni'ihau but much higher on Hawai'i. For coarse-scale rugosity, slope exerted the greatest influence on Hawai'i and the lowest on Ni'ihau and O'ahu (Figure 6b). The opposite was true for water depth. Here we note the apparent positive relationship between nearshore development and reef rugosity off Ni'ihau: Investigation of the mapping data indicated a highly localized spatial correlation between a small area of development and reef presence and rugosity relative to other sandy and cliff-like access points to the coastline.
Examination of model partial dependency indicated that water depth and reef slope had a positive linear relationship with finescale rugosity at the scale of the eight MHI (Figure 7). However, coarse-scale rugosity showed an asymptotic relationship with depth and slope (Appendix Figure 4). Moreover, while average fine-scale rugosity across all islands was fairly flat relative to distance to coast and all other factors (Figure 7), distance and wave power did show relationships with coarse-scale rugosity (Appendix Figure 4). We note an especially poor relationship between sea surface temperature and rugosity at either resolution.

DISCUSSION
Reef complexity, as indicated by maps of rugosity at 2 and 6 m resolutions, revealed a nested set of 3D habitat patterns across and within the eight MHIs. Coarse-scale rugosity was largely driven by reef slope and secondarily by depth; the opposite was true for fine-scale rugosity (Figure 6). Taken together, these two indicators suggest that geologic stage dominates the broader habitat structure of reefs across the archipelago. Younger islands harbor reefs on steep basaltic slopes, often with rapid drop-offs to extreme depths (Fletcher et al., 2008). During their formation via lava flows, these rock-dominated reefs resemble complex canyon, spire, tube, and spur-and-groove type patterns. As islands age, subside and undergo surficial erosional processes, nearshore shallow (<22 m depth) slopes become smoother (Fletcher et al., 2008), are reduced in rugosity, and show weakening relationships with coarse-scale rugosity (Figure 6).
In contrast to coarse-scale rugosity, the distribution and patterning of fine-scale features were far more driven by water depth than by reef slope (Figure 6), where inter-island differences in rugosity track differently with island age and stage. For example, the youngest island, Hawai'i, has the steepest nearshore slopes, where depth has relatively little effect on fine-scale rugosity as compared to O'ahu. For the latter, fine-scale rugosity tracks water depth, largely due to the prevalent role of live coral along forereef locations in the 15+ m range (see Asner et al., 2020b). Shallower inshore areas on O'ahu are largely devoid of rock or coral needed to drive fine-scale rugosity patterns. It is notable that spatial variation in fine-scale rugosity peaks in the middle-aged islands associated with the conglomerate geologic structure of Maui Nui (Figure 5), comprised of the islands of Maui, Lana'i, Kaho'olawe, and Moloka'i (Price and Elliott-Fisk, 2004). On these islands, reefs undergo large spatially explicit swings in fine-scale rugosity, which we hypothesize to be associated with highly variable coral habitat conditions moderated by substrate availability, currents, depth, and other sub-regional factors (Asner et al., 2020b). Such high-frequency variability in fine-scale rugosity is far less prevalent on older (Kaua'i, O'ahu) and younger (Hawai'i) islands (Figure 5).
Overlain on these geologically driven patterns, there may be an emerging set of 3D habitat conditions associated with nearshore development and other direct human impacts, but the signal remains variable by island (Figure 6). On islands with developed land-based infrastructure like O'ahu, fine-scale rugosity appears suppressed in areas of nearshore development (Figure 7), perhaps related to the suppression of live coral cover in these locations (Asner et al., 2020b). On other less developed islands like Hawai'i, fine-scale rugosity was positively related to nearshore development, perhaps because the relationship between land and reef patterns has not yet been sufficiently altered (Wedding et al., 2018). Repeat flights are required to identify areas with changing rugosity for subsequent mitigation actions related to nearshore development.
Beyond the drivers of reef complexity, our rugosity maps have revealed a nested set of patterns ranging in scale from local to archipelago. In a recent study across the MHI, the rugosity maps presented here were a primary determinant of resource fish, herbivore fish, and total fish biomass (Donovan et al., 2020). They also showed that the recovery potential for reef fish was strongly mediated by reef rugosity. These findings amplify both the importance of high-resolution reef complexity mapping and monitoring as well as the way that management can integrate habitat mapping into plans and actions to protect and restore reef fisheries (Nowlis and Friedlander, 2005;Wedding et al., 2008). FIGURE 7 | Partial dependence plots indicating the relationship between fine-scale rugosity and the top nine variables presented in Figure 6. Colored lines indicate relationships by island, and thick black lines indicate relationships for all islands combined. Coarse-scale rugosity is provided in Appendix Figure A4.
Two major limits to our study involve spatial resolution and maximum achieved depth. In our companion study describing the methods and validation for mapping 3D habitat complexity used here, Asner et al. (2020a) found that rugosity maps could be generated at 40-60 cm resolution, thereby revealing the location of smaller coral colonies and rock outcrops. These smaller features may be key to understanding a host of other reef patterns and processes such as fish and invertebrate dynamics, coral larval settlement, and disease (Hata et al., 2017;Eggertsen et al., 2020). The current study cannot resolve these ultrafine rugosity issues due to the practicalities of mapping an entire archipelago at high altitude and thus lower (2 m) spatial resolution.
Our maximum achieved mapping depth was 22 m (72 fsw). While much of the biological diversity of Hawaiian reefs is contained within our mappable depth range, critically important habitat in the lower euphotic and mesophotic zones are missed by our current approach (Rooney et al., 2010;Baldwin et al., 2018). While we continue to push the limits of spectroscopybased depth and rugosity mapping, with some indications of 30 m penetrability in clear waters (G.P. Asner, unpub. data), this is not likely to become operational in the near future due to the limits of solar photon flux at such depths paired with limits of our high-fidelity imaging spectrometers. While these current spectrometer-based approaches are limited to a very few manned aircraft platforms, the technology is downsizing and could become UAV (drone) based in the coming years. Additionally, these types of high-fidelity spectrometers are being prepared for low Earth orbit, such as the NASA Surface Biology and Geology mission (Schneider et al., 2019). Although the planned spatial resolution of 30 m will remain coarse for reef applications, the quality of the spectroscopic measurement will be analogous to the airborne capability today, perhaps allowing for coarse-scale rugosity mapping worldwide with the ecological fidelity proven here at an archipelago scale.
We found that water depth and reef slope are dominant drivers of both fine-and coarse-scale rugosity, where the greatest variation in rugosity occurred on middle-aged islands (e.g., Maui) in comparison to young and older islands (e.g., Hawai'i, Kaua'i). Nearshore development is an emerging secondary driver of decreasing rugosity, for example on O'ahu, the most highly developed island. Our rugosity maps are a tool that can be used by fisheries management and reef practitioners to prioritize areas for conservation and to design marine managed areas for greater resilience and recovery potential.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in Asner et al. (2020c).