Measuring Organization of Large Surficial Clasts in Heterogeneous Gravel Beach Sediments

<jats:p>The natural stratification and interlocking “organization” of armored sediments in heterogeneous, coarse-grained, beaches provides protection and enhances habitat for borrowing sedentary megafauna and macrofauna such as hard-shelled clams. Here, we develop a novel metric for quantifying sediment organization of large surficial beach clasts through sedimentologic and photogrammetric analyses of 37 lower intertidal heterogeneous gravel beaches in western Prince William Sound, Alaska (USA). Grain size, photogrammetric, and Wolman Pebble Count clast-size data from 64, ~1-m2 study plots are combined into a clast-size-independent “Organization Metric” to quantify the degree of organization in the meshed arrangement of larger surficial sediments. This metric was validated through field manipulation experiments and comparisons of adjacent plots characterized by different clast sizes. Application of this metric to subsets of Prince William Sound beaches that underwent differential treatment following the Exxon Valdez oil spill reveals persistent physical effects of artificial beach disturbance even 21 years after the cleanup. This has important implications for beach management (e.g., cleaning or dredging) and for the diverse and productive sedentary megafaunal assemblages that live within these sediments. Overall, this study provides a new approach for quantifying organization of heterogenous coarse sediments in diverse natural settings; in particular, heterogenous gravel beaches.</jats:p>

section_content: Heterogeneous coarse-grained beaches comprise those dominated by a mixture of gravel (intermediate-axis diameter of 2 mm or greater [1] ), and sand. They are generally composed of 15-68% sand (0. 063-2. 0 mm diameter) and 15-20% pebbles (4-640 mm diameter) or coarser material [2] [3] [4]. Although found throughout the world, they are most common along mountainous, tectonically active coasts and coasts proximal to presently or formerly ice-covered terrain (i. e., high-latitude paraglacial coasts [5] [6] [7]. For example, such gravelly and heterogeneous beaches comprise about one-third of the shoreline of the United Kingdom [8, 9] and many areas along the Gulf of Maine coast [10]. Coarse-grained beaches are also common along the paraglacial, leading-edge coast of southeast Alaska, where nearly 45% of total shoreline is either gravel or heterogeneous beaches [5]. 

A common feature of such heterogeneous coarse-grained beaches is the physical structure of surface gravel "armoring" [6, [11] [12] [13]. This attribute often has two components: (1) primary armoring, which characterizes the development of the two-layer sedimentary structure in which a coarse upper layer sits atop a finer, more heterogeneous sedimentary layer; and (2) organization, which reflects the stability of the armored surface, as determined by the orientation and/or interlocking nature of clasts in a relatively low-energy state (that is, with long axes closer to parallel to the beach slope). The latter can be characterized in terms of the orientation of the long (a), intermediate (b), and short (c) axes of a clast (see Table 1 ) with respect to the bed. Petrov [14] observed that the coefficient of flatness (c 2 /ab) for the pebbles and cobbles increased when coarse-bed material was exposed to wave action, demonstrating that surficial clasts became more "organized", i. e., armored clasts became reoriented to a lower, more stable, energy state, in response to reworking by waves. During the organization process, the upper-most coarse clasts can interlock (Figure 1 [15] ) thereby forming a loose "fabric," which sequesters fine sediments and organic matter in the underlying sedimentary matrix. Hayes and Michel [2, 15] ). In this case, intermediate-sized and some small clasts are removed from the upper bed of the intertidal beach by wave swash, whereas some small clasts are shielded by clasts which are too large to be moved. After initial armoring, further protection is driven by overpassing (intermediate-sized clasts bounce over armor and continue up-beach) and kinetic sieving (filtering of smaller particles into the interstices of larger particles) during both wave swash and backwash. Armoring may be a product of both models. 

Refined exploration of the ecological implications of sedimentological armoringand any natural or anthropogenic disturbance of that physical state-requires a quantitative means of describing and comparing the degree of clast organization on heterogenous gravel beaches. There has been substantial work seeking to quantify the stability of gravel surfaces, including the characterization of gravel bed roughness along riverine gravel (a) Segregation and armoring caused by a differentiation by clast shape and size during wave swash, followed by inverse grading due to progressively declining flow driven by infiltration of backwash (modified from Hayes and Michel [5]. (b) Armoring caused by storm wave swash flow over a mixed beach composed of clasts with different thresholds for movement (modified from Hayes and Michel [2, 15] ). In this case, intermediate-sized and some small clasts are removed from the upper bed of the intertidal beach by wave swash, whereas some small clasts are shielded by clasts which are too large to be moved. After initial armoring, further protection is driven by overpassing (intermediate-sized clasts bounce over armor and continue up-beach) and kinetic sieving (filtering of smaller particles into the interstices of larger particles) during both wave swash and backwash. Armoring may be a product of both models. 

Surface gravel armoring may play a critical role in coarse-grained ecosystems, particularly in lower intertidal and subtidal zones. It provides protection and substrate for diverse epibiotic assemblages of invertebrates and luxuriant macroalgae, which produce more organic material than almost any other intertidal habitat [16, 17]. It has also been suggested that gravel organization (the rotation and interlocking of clasts in the upper armor layer into a stable state) may be a prerequisite for development of perennial sedentary megafaunal assemblages living in those sediments, as well as for recovery following disturbance [18]. 

Refined exploration of the ecological implications of sedimentological armoring-and any natural or anthropogenic disturbance of that physical state-requires a quantitative means of describing and comparing the degree of clast organization on heterogenous gravel beaches. There has been substantial work seeking to quantify the stability of gravel surfaces, including the characterization of gravel bed roughness along riverine gravel point bars (e. g., Nikora et al. [19] ) and the armoring of gravel bed surfaces using flumes (e. g., Aberle and Nikora [20] ). However, much of these results are non-transferrable to coastal settings such as the heterogenous gravel beaches of Prince William Sound (PWS), because of the vastly different physical conditions. Fluvial gravel point-bar deposition and formation of surface fabrics are largely a product unidirectional high velocity, steady flow. Such fabrics are homogeneous, often composed of finer gravel with imbricated Gaussian particle distribution. In contrast, gravel beaches experience bi or multi-directional flow associated with the up-rush and backwash of breaking waves and longshore tidal currents. Additionally, the maximum physical conditions imparted to beach surface are constantly changing due to the rise and fall of the tides and varying wave conditions. Furthermore, beach sediments in PWS are composed of heterogeneous coarse lognormally distributed gravel and are not imbricated [5]. Moreover, most fabric-related gravel surface studies are aimed at understanding the hydrodynamic interrelationships between the gravel and stream flow conditions. In contrast, this study focused on how organization of the surface gravel imparts stability to the underlying sediment. Given these differences, we develop here a more computationally efficient approach, relying on surficial clast size and topographic elevation variability to describe the "organization" of the armoring clasts. This new clast-size-independent functional statistic (an "organization metric" or "OM") facilitates quantification of organization of large clasts concentrated at the surface of the heterogeneous gravel sediment. Though it may not be fully applicable to gravel riverbeds, it nonetheless allows for thorough and reproducible comparisons of the degree of armoring and organization in intertidal settings with widely varying sizes of armor clasts. 

We first present the conceptual framework of our approach for OM. We then detail development and field validation of this metric through sedimentologic and photogrammetric analysis of lower intertidal beaches in PWS. Finally, to demonstrate the utility of our new metric, we include an example of its use in exploring the alteration of clast organization within a subset of PWS beaches following the 1989 Exxon Valdez Oil Spill (EVOS).
section_name: Introduction

section_content: Prince William Sound is located along the northern margin of the Gulf of Alaska (Figure 2 ), a steep, bedrock-dominated coast, lying adjacent to the Alaska, Chugach, Talkeetna, and Wrangell-St. Elias mountain ranges. The protected and heavily dissected shoreline of PWS produces widely varying fetch-limited wave conditions. This characteristic, coupled with variable, rapidly fluctuating wind directions, reduces the ability of strong winds to generate high wave energy at some beaches. The result is a diverse beach setting with a wide range of hydrodynamic exposures [15]. Small pocket beaches that dominate PWS have longshore-varying textures and are composed of sediment largely derived from small streams, glacial deposits, and erosion of headlands composed of Tertiary volcanic and sedimentary rocks [21]. These beaches are relatively coarse and dominated by gravel and mud, with little intermediary sand. Hayes et al. attribute this characteristic to the predominance of shale and basalt, which break down to gravel-sized rock fragments and fine sediment [5]. Beach cusps, gravel imbrication, steep beach faces and cross-shore changes in clast shapes are all rare or absent on PWS beaches [5]. Except proximal to actively eroding cliffs, sediment textures generally coarsen in a seaward direction, across the macrotidal (tidal range: 4-5 m) beaches [22], with the lowest sections dominated by tectonically uplifted bedrock covered by a veneer of cobbles and boulders [15]. Armoring is common and tends to develop in the flatter middle and lower portions of PWS beaches [2, 5, 15, 23]. 

setting with a wide range of hydrodynamic exposures [15]. Small pocket beaches that dominate PWS have longshore-varying textures and are composed of sediment largely derived from small streams, glacial deposits, and erosion of headlands composed of Tertiary volcanic and sedimentary rocks [21]. These beaches are relatively coarse and dominated by gravel and mud, with little intermediary sand. Hayes et al. attribute this characteristic to the predominance of shale and basalt, which break down to gravel-sized rock fragments and fine sediment [5]. Beach cusps, gravel imbrication, steep beach faces and cross-shore changes in clast shapes are all rare or absent on PWS beaches [5]. Except proximal to actively eroding cliffs, sediment textures generally coarsen in a seaward direction, across the macrotidal (tidal range: 4-5 m) beaches [22], with the lowest sections dominated by tectonically uplifted bedrock covered by a veneer of cobbles and boulders [15]. Armoring is common and tends to develop in the flatter middle and lower portions of PWS beaches [2, 5, 15, 23]. All beaches in this region are geologically "young": most of the modern intertidal zone was shifted into its current vertical position through uplift of up to 11 m or several meters of downthrust associated with the 1964 Great Alaska Earthquake; many others were heavily disturbed by the ensuing tsunami [25]. Many of the gravel beaches of PWS were further artificially disarranged following the March 1989 EVOS. The discharge of 41 million liters of crude oil near the entrance to the Valdez Arm of PWS resulted in oiling of more than 780 km (~16%) of the intertidal zone in central and western PWS [26], an estimated 24% of which was characterized by heterogeneous gravel beaches [27]. A common postspill cleaning method employed was high-pressure/hot-water (HP/HW) washing, during which sediments were flushed with ~60 • C water applied at 50-120 psi and 950-1900 liters per minute. This often extended into the biologically productive lower intertidal zones, and commonly resulted in the removal of fine sediments (sand, silt, and clay) from shallow beach sediments and disrupted the orientation of large clasts on these beaches (DCL, pers. obs. ; Figure 3 ). 

post-spill cleaning method employed was high-pressure/hot-water (HP/HW) washing, during which sediments were flushed with ~60 °C water applied at 50-120 psi and 950-1900 liters per minute. This often extended into the biologically productive lower intertidal zones, and commonly resulted in the removal of fine sediments (sand, silt, and clay) from shallow beach sediments and disrupted the orientation of large clasts on these beaches (DCL, pers. obs. ; Figure 3 ). We selected 37 beaches from across western PWS, which represented a combination of those which had, and had not, been subjected to disturbance associated with HP/HW washing following the EVOS, as per ADNR [28] (Figure 2 ). Most are located in protected embayments and small coves where sediments are primarily composed of a mixture of small boulders, cobbles, pebbles, sand, and silt (i. e., heterogeneous gravel habitats). However, several relatively more exposed beaches such as LA16 and LA18 on Latouche Island and the northern sites on Disk Island were dominated by large cobbles and small boulders. Longshore textural diversity within our beach sites was common, with some beaches dominated in some areas by pebbles to cobbles and in others by cobbles to boulders. Thus, most sites were subdivided into two smaller study areas, chosen to reflect this textural diversity. Sampling was conducted within the lower intertidal zone, at elevations ranging from -0. 6 m to +0. 9 m Mean Lower Low Water. We selected 37 beaches from across western PWS, which represented a combination of those which had, and had not, been subjected to disturbance associated with HP/HW washing following the EVOS, as per ADNR [28] (Figure 2 ). Most are located in protected embayments and small coves where sediments are primarily composed of a mixture of small boulders, cobbles, pebbles, sand, and silt (i. e., heterogeneous gravel habitats). However, several relatively more exposed beaches such as LA16 and LA18 on Latouche Island and the northern sites on Disk Island were dominated by large cobbles and small boulders. Longshore textural diversity within our beach sites was common, with some beaches dominated in some areas by pebbles to cobbles and in others by cobbles to boulders. Thus, most sites were subdivided into two smaller study areas, chosen to reflect this textural diversity. Sampling was conducted within the lower intertidal zone, at elevations ranging from -0. 6 m to +0. 9 m Mean Lower Low Water.
section_name: Study Area

section_content: We developed, tested, and validated our new organization metric through analysis of sediments in the lower intertidal zone of 37 heterogenous, coarse-grained beaches in western PWS. Specific methods employed include: (1) textural analyses of sediment samples collected along shallow depth profiles at each beach; (2) photogrammetric analyses of 1-2 three-dimensional plots (~1 m 2 ) within the lower intertidal zone at each of these beaches (64 total photoplots); and (3) Wolman Pebble Counts [29] of larger clasts in the surficial fabric of these plots.
section_name: Materials and Methods

section_content: Surficial gravel armoring and organization can occur in response to both uni and bidirectional currents [19]. Early studies focused on the former, quantifying armoring of mixed gravel beds in laboratory flumes or mountain streams using roughness terms [30], an approach still widely used today [31]. For example, Nikora et al. [19] measured gravel-bed roughness as a random field of vertical measurements (z values) in a three-dimensional space. They determined "roughness" by evaluating the variation of the population of z values (i. e., their standard deviation, σ z ) in a relatively uniform gravel bed, and evaluated these against the axial dimensions (i. e., a, b, and c axes) of the sampled clasts. However, this approach fails to characterize the three-dimensional nature of surface organization in sediments with variable textures, in which processes such as kinetic sieving produce a surface layer that is generally coarser than underlying deposits [32]. Aberle and Nikora [20] used random field elevation based on laser scanning techniques in laboratory flumes to develop Digital Terrain Models (DTMs) of bed surfaces that describe the structure or roughness of organized gravel sediments on scales ranging from several millimeters to that of the largest clasts within the bed. They found that vertical roughness of gravel beds can be defined by the standard deviations of bed elevations, and that these bed elevations can be used to accurately quantify the development of a coarse surface layer of clasts within that layer. However, this approach is difficult to translate to the field, and does not allow for quantification and comparison of surface clast orientations across a range of clast sizes. 

Surficial sediment organization represents the degree to which surface clasts are oriented in their most kinetically stable state. At maximum organization, the long (a) and wide (b) axes of large clasts are parallel to the surface of the bed and the shortest (c) axis is oriented orthogonally to the bed slope. The farther clasts protrude from the surface of the beach matrix, the more likely that those clasts are out of their most stable position; that is, the less organized the surficial sediment. This measure could be characterized by the range of "elevations" above the matrix bed level (z values), quantified as the mean z value (zbar) or-following Nikora et al. [19] -the surface roughness (σ z ) of a given bed (Table 1 ). It follows that more highly organized surface sediments would, on average, have lower values of zbar and σ z ; that is, these would have a lower average elevation above the sedimentary matrix level and less elevation variability (Figure 4 ). The maximum value of zbar for a sediment plot would be approximately equal to the mean long-axis dimension of clasts within that plot (abar; Table 1 ) if (1) all large clasts were of similar size (same length of long ['a'] axes), (2) those large clasts were all positioned with their long axes oriented perpendicularly to the beach slope, and (3) the base of the largest clasts were at the zero-elevation plane. This would represent the minimum theoretical degree of surface clast organization. At sites exhibiting progressively more organization, an increasingly greater percentage of clasts would have their long axis oriented closer to horizontal, and thus be in a more "stable" state. 

However, this conceptual framing implies that all clasts are of similar or identical size, both within a plot and across the population of sites being compared. In heterogenous sediments, the range of surface elevations (z values) within a plot may reflect, in part, the range in clast sizes, irrespective of the orientation of those clasts. For example: for two beaches in which all clasts are aligned in their lowest-energy configuration (a and b axes aligned with the plane of the beach slope), the surface rugosity (zbar, or σ z ) would be greater at the beach with coarser clasts (Figure 4, second panel). This reflects that the c-axis of those coarser clasts is larger-and thus the elevation difference between surface peaks and valleys would be greater-than in a plot with similarly oriented, but overall smaller, clasts (Figure 4, bottom panel). Thus, measures of roughness would depend, in part, on clast size. Moreover, given similar input energy (for example, from breaking waves on a beach), smaller clasts would be expected to be more easily or quickly rotated into a lower-energy, bed-parallel orientation. Thus, within a shoreline reach with consistent hydrodynamic energy, beach sub-sections characterized by finer clasts should become more organized-or organize more rapidly-than those consisting of coarser clasts. 

This conceptualization highlights the importance of accounting for clast size in characterizing surficial organization, thereby allowing for meaningful inter-site comparisons regardless of bed clast size. In his analysis of streambeds, Wolman [29] used measurements of the b-axis of individual clasts to characterize the diversity in surficial clast sizes. This axis provides a measurement for large clasts that is analogous to mean grain size (Mz) for the bulk sediment as measured through sieve analysis. In contrast, the level of organization of clasts in a bed is predominantly a function of the orientation of the longest (a) axis of the large clasts: the longer the clasts, the greater maximum surface rugosity in a fully disorganized bed (Figure 4, top panel). We thus consider the average length of the long-axis of clasts in a given plot (abar) to provide the most suitable measure by which to account for clast-size variability. 

Following this theoretical framework, we designed a study to quantitatively refine and validate a new organization metric (OM), relying on field measurements of both surface roughness and clast size. In the following section we present this approach, the resulting OM, and validation tests for this new metric. However, this conceptual framing implies that all clasts are of similar or identical size, both within a plot and across the population of sites being compared. In heterogenous sediments, the range of surface elevations (z values) within a plot may reflect, in part, the range in clast sizes, irrespective of the orientation of those clasts. For example: for two beaches in which all clasts are aligned in their lowest-energy configuration (a and b axes aligned with the plane of the beach slope), the surface rugosity (zbar, or σz) would be greater at the beach with coarser clasts (Figure 4, second panel). This reflects that the caxis of those coarser clasts is larger-and thus the elevation difference between surface peaks and valleys would be greater-than in a plot with similarly oriented, but overall smaller, clasts (Figure 4, bottom panel). Thus, measures of roughness would depend, in part, on clast size. Moreover, given similar input energy (for example, from breaking waves on a beach), smaller clasts would be expected to be more easily or quickly rotated into a lower-energy, bed-parallel orientation. Thus, within a shoreline reach with consistent hydrodynamic energy, beach sub-sections characterized by finer clasts should be-
section_name: Conceptual Framework for Quantifying Organization in Heterogenous Gravel Beaches

section_content: At each site, bulk sediment samples were collected from the side of shallow (<15 cm deep) infaunal excavation by inserting a 15-cm long, 10-cm diameter core sampler horizontally into the undisturbed sediment of the excavation wall. Samples collected in this manner provide a representative profile of grain size in three 5-cm deep horizons, i. e., upper, middle, and lower strata. Care was taken to include all coarse clasts that could fit into the core sampler. Where present, clasts with b-axes ≥10 cm were excluded due to size limitations of the core sampler. 

Coarse fractions (>4 phi [ϕ] [0. 0625 mm]) were dry sieved at 1-ϕ intervals following procedures of Folk and Ward [33]. Fine fractions separated by wet sieving were analyzed using a Beckman Coulter LS 13320 Laser Diffraction Particle Size Analyzer to 11ϕ (0. 00049 mm; fine clay). Grain-size statistics were calculated following Folk and Ward [33] using the logarithmic (inclusive) graphical method.
section_name: Granulometric Analyses

section_content: General morphological characteristics, including overall beach slope, sediment composition, and fetch exposure were qualitatively described for each study site. At each, plots, roughly square and ~1 m 2 in area, were marked out using a series of 10-cm diameter coded targets (Figure 5 ). The surface of each plot was cleared of all plant cover and large motile invertebrates. Care was taken to ensure that clasts within these plots remained undisturbed during the clearing process. Vegetation was cut from the individual clasts, exposing the surface of the clasts such that morphological analysis captured the shape and form of the bare clasts (Figure 5b, c ). Lateral reference scales within each plot were provided by a coded meter stick placed on the beach surface (Figure 5b ). Vertical reference was provided by a PVC pipe oriented with a hand level (Figure 5a ). A series of high-resolution stereophotograph pairs was taken of each plot (e. g., Figure 5 ) using a Nikon D90 camera with a 20. 0-mm lens specially calibrated for working with the PhotoModeler Scanner TM (EOS ® Systems Inc., Vancouver, BC) software package. Paired photos were offset laterally by <50 cm. These pairs were taken from at least two angles and elevations from each side of a given plot to ensure all visible surfaces of each clast were photographed. This provided stereophoto pairs from a variety of orientations (Figure 5b, c ). A series of 5-10 stereo pairs was collected from each plot. Generally, plots dominated by coarser-grained gravel required a greater number of photos than those with finer clasts due to the number of orientations needed to capture the three-dimensional morphology of larger cobbles and boulders. 

Following the collection of the photographs, Wolman Pebble Counts [32] were conducted for each plot in order to determine the long (a), intermediate (b), and short (c) axes for the largest 15-47 surface clasts within each plot. The range in the number of clasts measured reflects the diversity in clast size within a given plot. Those with larger and more diverse clast sizes required additional clast counts to characterize clast diversity within the plot. These were averaged for each plot and are given as abar, bbar, and cbar, respectively. 

Finally, in an experiment designed as a validation test for the method and the organization metric, the surface of four randomly selected plots was manually disturbed following completion of pebble counts. Clasts comprising the surface armor were vigorously disturbed with an entrenching tool, shovel, or trowel. The photo plots were then re-photographed and analyzed in the same manner as for the initial, undisturbed plots.
section_name: Field Geomorphology and Photogrammetry

section_content: Stereo-photo pairs were analyzed using PhotoModeler Scanner TM. This system produces a dense point cloud from photographs of textured surfaces (surface sediment plots, in this case). Photos were inputted into the software package and stereopairs identified. Coded targets, visible on all plot photos, were manually identified on each photo and cross-referenced between photos. Photos were then "idealized" to account for artifacts associated with the curvature of the camera lens. The resulting photos were edge-warped such that each photo captured a true projection of space, without lens distortion and with square pixels and a centered principal point. Photo scale and orientations were calibrated using the meter stick and the vertical pipe. 

Three-dimensional maps of plots were created at ~1-mm horizontal resolution (Figure 5d ) and trimmed to the given plot areas using coded targets as plot area boundaries. Each point within a plot is assigned an x, y coordinate in non-referenced space, and an elevation value (z) above an artificial horizontal plane. The number of coordinate points obtained from a given 1-m 2 photoplot ranged from ~250,000 to ~700,000. This variability reflects differences in the area covered by the photoplots and the complexity of the surficial fabric formed by the clasts. The resulting point mesh data were exported as digital terrain models (DTMs) into MATLAB for comparative analyses. All data underwent first-order detrending to remove the slope within the data set that resulted from the natural slope of the beach. The resulting z values were corrected such that the minimum elevation point in a given plot was assigned a value of 0 mm. Surface elevation (z-axis) statistics such as maximum, minimum, mean (zbar), median, and standard deviation (σ z ) elevations were calculated for each plot. Horizontal slices were then taken through the DTMs at 5-mm increments, at the median elevation, and at elevations representing the 10th and 25th percentiles above the minimum elevation in a given plot. The percent of data points above and below each of those slices was then calculated. 

Combined with our Wolman Pebble Counts and surface sediment grain-size analysis, these data provide variables that were used to characterize size, and surficial organization of sediments within each plot.
section_name: Processing and Analysis of Photogrammetric Data

section_content: Site comparisons used for validation and testing of our approach were evaluated using resampling techniques for comparison of means, analysis of variance (ANOVA), and regression analyses. These were calculated using a Microsoft Excel add-in developed by Resampling Stats [35]. Two-factor ANOVA was calculated using the StatPlus add-in to Excel. Where particular differences were expected based on field observations, 1-tailed tests were used. 

For comparing regression slopes, the paired linear regression approach outlined in Zar [36], was implemented with IGOR Pro version 6. 22A (http://www. wavemetrics. com; accessed on 9 March 2022). This statistic compares both the slope of two regressions and the mean y-values of the data comprising two regressions. First, it assesses whether rates of change for the slopes of two regressions are the same. Next, it compares the "statistical elevations" of two populations (as distinct from z-values in the photogrammetric analyses) to assess whether they are the same. A significant difference in either comparison indicates that the populations of values are statistically different. 

Finally, given the heterogeneity of PWS lower intertidal beaches in terms of lithology and wave exposure (particularly in comparison to conditions experienced in laboratory settings), and to minimize false negatives, we adopted an α value of 0. 1 as the critical level of significance for all statistical testing. Values of "p" between 0. 1 and 0. 2 are considered to represent trends. Standard deviation (SD) is used to represented degree of variability for individual measurements (e. g., OM) and standard error (SE) is used where variables are means (e. g., zbar or abar).
section_name: Inferential Statistics

section_name: Results

section_content: Treatment status, nominal tidal elevation, mean clast length (abar) calculated from Wolman Pebble Counts, and photoplot data (elevation standard deviation (σ z ), and the calculated values for our organization metric (OM)) are presented in Table 2. Mean particle size (Mz) and % cobble/pebble/granule (gravel) from grain-size analyses of bulk-sediment samples from the upper sediment horizon at each site, as well as its nominal tidal elevation, are presented in Supplemental Table S1. Sediment textures varied considerably among sites. The mean grain size (Mz) of the upper 5 cm of sediment for all 37 sites surveyed ranged from 3. 01 to 51. 85 mm and averaged 14. 7 ± 10. 3 mm (Supplemental Table S1 ). These values reflect the presence of a very wide range of particle size classes in the sediments among the sites, i. e., clay to small cobble. Percent gravel averaged 75. 3% ± 16. 3 and ranged from 32. 1 to 99. 9%. 

Among surficial clasts in each plot, the subset of the largest clasts measured through Wolman Pebble Counts had a mean long-axis dimension (abar) of 4. 8 to 19. 3 cm, averaging 9. 19 ± 0. 41 cm across all sites (Table 2 ). Calculated abar and Mz values were not correlated (p > 0. 5; resampling linear regression), indicating that the length of surficial clasts was not related to the composition of the underlying whole sediments. 

We found that both Mz and percent gravel were highest in the upper stratigraphic horizon and were lower in the middle and lower horizons (Table 3 ). Sand content was highest at the middle horizon. Silt/clay content was lowest at the upper level and highest in the lower horizon. The differences between the horizons were significant for all variables except the percentage of silt/clay. The mean surface elevation of sediments in our slope-corrected photoplots (zbar) was, on average, 4. 51 ± 0. 25 cm above the artificial horizontal plane at the lowest imaged elevation (z = 0; see Figure 2 ). This value is ~50% of mean measured abar for these sites. The standard deviation of these photoplot elevations (σ z ) averaged 1. 51 ± 0. 11 cm and ranged from 0. 52 to 6. 09 cm. As expected, both zbar and σ z exhibited very strong correlations with abar values (Figure 6 ): The correlation of σ z to abar (p << 0. 0001; resampling linear regression) was stronger than that for zbar (p = 0. 0001) and thus provides a stronger statistic with which to evaluate clast organization. However, σ z is not clast-size independent, as it exhibits a very strong relationship with abar (Figure 6 ). Thus, as anticipated, σ z is not suitable by itself for describing organization among sites with varying clast sizes. 

highest at the middle horizon. Silt/clay content was lowest at the upper level and highest in the lower horizon. The differences between the horizons were significant for all variables except the percentage of silt/clay. 3. 0 ± 3. 1 4. 0 ± 6. 0 6. 7 ± 14. 9 0. 191

The mean surface elevation of sediments in our slope-corrected photoplots (zbar) was, on average, 4. 51 ± 0. 25 cm above the artificial horizontal plane at the lowest imaged elevation (z = 0; see Figure 2 ). This value is ~50% of mean measured abar for these sites. The standard deviation of these photoplot elevations (σz) averaged 1. 51 ± 0. 11 cm and ranged from 0. 52 to 6. 09 cm. As expected, both zbar and σz exhibited very strong correlations with abar values (Figure 6 ): The correlation of σz to abar (p << 0. 0001; resampling linear regression) was stronger than that for zbar (p = 0. 0001) and thus provides a stronger statistic with which to evaluate clast organization. However, σz is not clast-size independent, as it exhibits a very strong relationship with abar (Figure 6 ). Thus, as anticipated, σz is not suitable by itself for describing organization among sites with varying clast sizes.
section_name: Field and Laboratory Results

section_content: Given the myriad clast sizes observed on PWS beaches, and the range in the degree of organization of those clasts in the lower intertidal zone, variability in elevations (z) within a plot can be large. We found a significantly stronger correlation between clast size (abar) and the range of elevations present in a given plot (i. e., σz; where r = 0. 67, p = <<
section_name: Derivation of the Organization Metric

section_content: Given the myriad clast sizes observed on PWS beaches, and the range in the degree of organization of those clasts in the lower intertidal zone, variability in elevations (z) within a plot can be large. We found a significantly stronger correlation between clast size (abar) and the range of elevations present in a given plot (i. e., σ z ; where r = 0. 67, p = << 0. 0001; Figure 6 ) than between abar and the average plot elevation (zbar; where r = 0. 48, p = 0. 0001; resampling regression). This analysis confirms our theoretical framework for applying a measure of surface roughness (as characterized by σ z ) in our organization metric and normalizing by the clast-size measure to which organization is most sensitive (abar). 

Following these considerations, we derive the following formula for calculating our dimensionless organization metric (OM):

in which both abar and σ z values are given in cm, so that higher values represent higher degrees of surficial organization (clasts situated in most stable orientations). 

The calculated OM values for all sites ranged from 2. 27 to 14. 06 (Table 2 ; Figure 6 ), with a mean value of 6. 89 ± 2. 307. The OM values are independent of abar (p = 0. 25). Thus, equal values of OM should indicate that the mean orientation (flatness) of the clasts relative to beach angle is similar at sites at which clast lengths differ substantially.
section_name: Derivation of the Organization Metric

section_content: We employed three independent tests to evaluate the OM function. Together, these demonstrate that differences in the OM values accurately reflect differences in surficial clast organization. 

The first test relies on the fact that, all other factors being equal, the distribution of z-values (the "elevation" at each x,y point in the DTMs) is itself an indicator of the degree of organization, as demonstrated by Aberle and Nikora [20]. In Figure 7, we present the mean percentage of z-values extending varying distances above the standardized z = 0 matrix bed elevation for a given DTM for all undisturbed plots and four disturbed plots surveyed (KN5-1, KN5-2, KN208-1, and KN208-2). Thus, at the 50-mm slice, 35% of the z-values for the undisturbed plots extend above the standardized zero elevation whereas nearly 73% of the z values for the four disturbed plots extend above that elevation, indicating a greater degree of disorientation in the clasts in the disturbed plots. This analysis demonstrates that a larger proportion of the z-values forming the surfaces of the DTM for the undisturbed plots occurred significantly closer to the standardized base level (zero on the x-axis) than in the disturbed plots (p << 0. 0001; resampling comparison of means). 

The second test involved an assumption that smaller clasts are more easily rotated into a flattened position because less hydrodynamic (wave or flow) energy is required to move them. Thus, assuming that seasonally averaged wave energy along a given short stretch of beach distal from headlands (<50 m) is approximately equal, areas of that beach with finer clasts (smaller abar) should be more organized than those characterized by coarser clasts (larger abar). It follows that OM values for the plot with finer clasts should, on average, be higher (more apparently well organized) than the one with coarser clasts. The predominance of this pattern is clear (Figure 8 ). At these sites, plots with smaller clasts had higher values of OM (were more organized) at 21 of 29 sites. Mean OM for plots with smaller clasts (7. 3 ± 0. 4) was significantly greater than that observed for plots with larger clasts (6. 3 ± 0. 5; Table 4 ; p = 0. 080; paired 1-tailed resampling comparison of means). 

Our third field functionality test for OM compared pre and post-disturbance values for OM and σ z for four plots that were manually disturbed with a shovel following the initial photo-documentation. We predicted that OM in the artificially disturbed plots would be lower than in those same plots prior to manual disturbance. Indeed, following disturbance, zbar increased 79% ± 51% (p < 0. 0001; paired 1-tailed resampling comparison of means; mean σ z also increased (52. 0% ± 40. 0%; p < 0. 0001; Table 5 ) and mean OM decreased significantly (31. 1% ± 8. 7%; p < 0. 0001; Figure 9 ). All three variables responded significantly to disturbance. 

of organization, as demonstrated by Aberle and Nikora [20]. In Figure 7, we present the mean percentage of z-values extending varying distances above the standardized z = 0 matrix bed elevation for a given DTM for all undisturbed plots and four disturbed plots surveyed (KN5-1, KN5-2, KN208-1, and KN208-2). Thus, at the 50-mm slice, 35% of the zvalues for the undisturbed plots extend above the standardized zero elevation whereas nearly 73% of the z values for the four disturbed plots extend above that elevation, indicating a greater degree of disorientation in the clasts in the disturbed plots. This analysis demonstrates that a larger proportion of the z-values forming the surfaces of the DTM for the undisturbed plots occurred significantly closer to the standardized base level (zero on the x-axis) than in the disturbed plots (p << 0. 0001; resampling comparison of means). The second test involved an assumption that smaller clasts are more easily rotated into a flattened position because less hydrodynamic (wave or flow) energy is required to move them. Thus, assuming that seasonally averaged wave energy along a given short stretch of beach distal from headlands (<50 m) is approximately equal, areas of that beach with finer clasts (smaller abar) should be more organized than those characterized by coarser clasts (larger abar). It follows that OM values for the plot with finer clasts should, on average, be higher (more apparently well organized) than the one with coarser clasts. The predominance of this pattern is clear (Figure 8 ). At these sites, plots with smaller clasts had higher values of OM (were more organized) at 21 of 29 sites. Mean OM for plots with smaller clasts (7. 3 ± 0. 4) was significantly greater than that observed for plots with larger clasts (6. 3 ± 0. 5; Table 4 ; p = 0. 080; paired 1-tailed resampling comparison of means). decreased significantly (31. 1% ± 8. 7%; p < 0. 0001; Figure 9 ). All three variables responded significantly to disturbance.
section_name: Validation of the Organization Metric

section_content: The conceptual framework for OM suggests that the mean orientation of the long (a) axis of armoring clasts with the plane of the beach (the "clast angle") should decline as OM increases, i. e., the clasts become more parallel with the surface of the sediment. We calculated the average clast orientation in each photoplot by designating abar as the adjacent side of a right triangle and zbar as the opposite side. The tangent of the ratio of these sides (zbar/abar), was used to approximate the average angle formed by the armoring clasts with the surface of the beach. The resulting comparison between OM and mean clast angle, shown in Figure 10, indicates that, indeed, mean photoplot clast angle decreases at higher values of OM, and that the relationship is highly significant (p<<0. 0001; resampling regression).
section_name: Relationship between OM and Clast Angle

section_content: The question arises: "How can an organization metric be used?" To demonstrate, we present the outcomes of a case study using our new OM to explore the long-term (21 years) effects of post-spill cleanup efforts on PWS beaches. Specifically, we compared the degree of surficial clast organization in the lower intertidal zone of beaches which were left untouched following the EVOS with those that likely experienced significant physical disturbance caused by HP/HW washing (designations of our sites shown in Figure 3 and given in Table 2 ). In addition to demonstrating the efficacy of our approach and OM, a goal of this case study is to assess whether those beaches which underwent HP/HW washing following the EVOS would be less organized than those of unwashed beaches. The likelihood of differing degrees of clast organization associated with post-spill treatment was first proposed by Lees and Driskell [18], who concluded that large sedentary megafauna, especially adult and sub-adult hard-shelled clams, were significantly less abundant in washed beaches. They suggested that reduced organization would affect the rate of recovery of assemblages that were rich in the areas prior to the EVOS cleanup. Comparing beach sedimentology below the surface through our depth-profile samples, we were further able to test washing effects extended deeper into the sediment than just disturbance of surficial armor. 

Mearns reports that only 89% of the PWS shoreline segments oiled during the EVOS were treated in some manner [37]. Of these, he reported that a large fraction (ca. 79%) was treated to remove oil with approaches such as HP-HW washing, which are very disruptive to surficial stratigraphy and organization (e. g., Figure 4 ). Unfortunately, records from the post-spill clean-up period are incomplete and do not allow for identification of all beaches or beach segments that were treated with such washing. As such, we specifically selected beaches and beach segments that had been personally observed to have undergone this treatment during the clean-up period ("observed washed"; DCL) or were highly likely to have undergone this treatment ("presumed washed") based on interrelationships between location, degree of oiling, and governmental agency records [28, 38]. We group these as "presumed" washed sites. For analyses in the following section, we have grouped the latter two categories and refer to them as O/P washed sites. 

Given uncertainty regarding post-EVOS treatments, we applied a 'weight-of-evidence' argument to the outcomes of our analyses for those sites classified as "O/P washed", in order to improve confidence in this designation. One can expect HP/HW washing to alter a variety of sedimentary and stratification characteristics, beyond the hypothesized impact on surficial organization. For example, for beaches that underwent HP/HW washing, we would expect to see predictable relationships, on average (Table 6 ):
section_name: Case Study: Use of OM to Quantify the Long-Term Impact of Beach Washing

section_content: The question arises: "How can an organization metric be used?" To demonstrate, we present the outcomes of a case study using our new OM to explore the long-term (21 years) effects of post-spill cleanup efforts on PWS beaches. Specifically, we compared the degree of surficial clast organization in the lower intertidal zone of beaches which were left untouched following the EVOS with those that likely experienced significant physical disturbance caused by HP/HW washing (designations of our sites shown in Figure 3 and given in Table 2 ). In addition to demonstrating the efficacy of our approach and OM, a goal of this case study is to assess whether those beaches which underwent HP/HW washing following the EVOS would be less organized than those of unwashed beaches. The likelihood of differing degrees of clast organization associated with postspill treatment was first proposed by Lees and Driskell [18], who concluded that large sedentary megafauna, especially adult and sub-adult hard-shelled clams, were significantly less abundant in washed beaches. They suggested that reduced organization would affect the rate of recovery of assemblages that were rich in the areas prior to the EVOS cleanup. Comparing beach sedimentology below the surface through our depth-profile samples, we were further able to test washing effects extended deeper into the sediment than just disturbance of surficial armor. 

Mearns reports that only 89% of the PWS shoreline segments oiled during the EVOS were treated in some manner [37]. Of these, he reported that a large fraction (ca. 79%) was treated to remove oil with approaches such as HP-HW washing, which are very disruptive to surficial stratigraphy and organization (e. g., Figure 4 ). Unfortunately, records from the post-spill clean-up period are incomplete and do not allow for identification of all beaches or beach segments that were treated with such washing. As such, we specifically selected beaches and beach segments that had been personally observed to have undergone this treatment during the clean-up period ("observed washed"; DCL) or were highly likely to have undergone this treatment ("presumed washed") based on interrelationships between location, degree of oiling, and governmental agency records [28, 38]. We group these as "presumed" washed sites. For analyses in the following section, we have grouped the latter two categories and refer to them as O/P washed sites. 

Given uncertainty regarding post-EVOS treatments, we applied a 'weight-of-evidence' argument to the outcomes of our analyses for those sites classified as "O/P washed", in order to improve confidence in this designation. One can expect HP/HW washing to alter a variety of sedimentary and stratification characteristics, beyond the hypothesized impact on surficial organization. For example, for beaches that underwent HP/HW washing, we would expect to see predictable relationships, on average (Table 6 ):

1. 

Mz in surficial sediments at unwashed sites should be finer than at O/P washed sites because flushing associated with HP/HW washing would have removed much of the fine fraction.
section_name: Case Study: Use of OM to Quantify the Long-Term Impact of Beach Washing

section_content: The proportion of gravel/cobble therefore should be greater, and the proportion of finer sediments would be lesser, than at unwashed sites.
section_name: 2.

section_content: Mz in the surficial 5-cm horizon should be coarser than in deeper horizons due to that same flushing of fine sediments; and 4. 

The proportion of silt/clay should be lower at all horizons at O/P washed sites than at unwashed sites. 

We would also expect to see predictable relationships in abar and the photogrammetric values, on average (Table 7 ):
section_name: 3.

section_content: Because washing should not affect the large clasts, abar should be similar.
section_name: 5.

section_content: zbar and σ z should be greater than at unwashed sites. 7. 

OM should be lower than at unwashed sites. Weight-of-evidence in these tables argues strongly that, on average, the two treatment categories are distinctly different, i. e., one group of sites was unwashed, and the other group was washed. Mz, which responds to loss of sand and silt/clay due to washing, differed significantly between unwashed (9. 00 ± 0. 95 mm) and O/P washed sites (18. 21 ± 2. 40 mm; p-0. 0076, resampling comparison of means; Table 6 ). Gravel was greater at O/P washed sites (83. 0% ± 2. 4 vs. 73. 9 ± 2. 3; p = 0. 004) because a large volume of sand and silt/clay that was washed away. Moreover, Mz in the lower horizon was significantly lower than in the upper horizons (p-0. 014; 2-tailed resampling ANOVA). Finally, % silt/clay was significantly lower in upper and lower horizons than at O/P washed sites (p << 0. 0001 and p = 0. 08, respectively) because it had been flushed out by washing. 

Notably, differences in abar between unwashed and O/P washed sites were not significant (p = 0. 94; Table 7 ). This was predictable because large clasts likely would not be added or removed from an area by washing. It was unexpected that zbar was not significantly higher at O/P washed sites (p = 0. 22). However, σ z was significantly higher at O/P washed sites (p = 0. 073) and OM was significantly higher sites at unwashed sites, as predicted above (p = 0. 06). 

Comparison of OM values for unwashed and O/P washed sites (Figure 11 ) reveals overall substantially lower degree of organization and considerably more variability in OM in the O/P washed sites. Slopes and separation on the y-axes (statistical elevations) of the OM regressions for the unwashed and O/P washed sites were compared using the paired linear regression approach (Figure 11 ). The data clouds overlap considerably, reflecting the great variability in environmental conditions at these sites. Nevertheless, the slopes are not statistically different (p > 0. 5; paired linear regression as calculated by IGOR), again demonstrating that OM fits a primary criterion for metric suitability (clastsize independence). However, despite the high variability among the sites, the statistical elevations are statistically different (p = 0. 025), with the OM regression for unwashed sites significantly higher than for O/P washed sites. However, both are relatively weak, as indicated by the relatively low p-values for the regressions (p > 0. 63 and p > 0. 34, respectively). This indicates that beach washing reduced armor clast organization at O/P washed sites and this effect persisted 21 years after the cleanup. 

Weight-of-evidence in these tables argues strongly that, on average, the two treatment categories are distinctly different, i. e., one group of sites was unwashed, and the other group was washed. Mz, which responds to loss of sand and silt/clay due to washing, differed significantly between unwashed (9. 00 ± 0. 95 mm) and O/P washed sites (18. 21 ± 2. 40 mm; p-0. 0076, resampling comparison of means; Table 6 ). Gravel was greater at O/P washed sites (83. 0% ± 2. 4 vs. /73. 9 ± 2. 3; p = 0. 004) because a large volume of sand and silt/clay that was washed away. Moreover, Mz in the lower horizon was significantly lower than in the upper horizons (p-0. 014; 2-tailed resampling ANOVA). Finally, % silt/clay was significantly lower in upper and lower horizons than at O/P washed sites (p << 0. 0001 and p = 0. 08, respectively) because it had been flushed out by washing. 

Notably, differences in abar between unwashed and O/P washed sites were not significant (p = 0. 94; Table 7 ). This was predictable because large clasts likely would not be added or removed from an area by washing. It was unexpected that zbar was not significantly higher at O/P washed sites (p = 0. 22). However, σz was significantly higher at O/P washed sites (p = 0. 073) and OM was significantly higher sites at unwashed sites, as predicted above (p = 0. 06). 

Comparison of OM values for unwashed and O/P washed sites (Figure 11 ) reveals overall substantially lower degree of organization and considerably more variability in OM in the O/P washed sites. Slopes and separation on the y-axes (statistical elevations) of the OM regressions for the unwashed and O/P washed sites were compared using the paired linear regression approach (Figure 11 ). The data clouds overlap considerably, reflecting the great variability in environmental conditions at these sites. Nevertheless, the slopes are not statistically different (p > 0. 5; paired linear regression as calculated by IGOR), again demonstrating that OM fits a primary criterion for metric suitability (grainsize independence). However, despite the high variability among the sites, the statistical elevations are statistically different (p = 0. 025), with the OM regression for unwashed sites significantly higher than for O/P washed sites. However, both are relatively weak, as indicated by the relatively low p-values for the regressions (p > 0. 63 and p > 0. 34, respectively). This indicates that beach washing reduced armor clast organization at O/P washed sites and this effect persisted 21 years after the cleanup. In Figure 12, we present the mean percentage of z-values extending varying distances above the standardized z = 0 (see Figure 2 ) elevation for a given sets of DTMs for all unwashed, O/P washed surveyed, and four plots manually disturbed during the study (KN5-1, KN5-2, KN208-1, and KN208-2). At the 50-mm slice, 31% of the z-values for the unwashed plots extend above the standardized zero elevation, compared to nearly 37% of the z-values for O/P washed plots. Finally, 73% of the z-values for the four disturbed plots extend above that elevation. This pattern continues for the remainder of the O/P washed and disturbed curves. This analysis demonstrates that a larger proportion of the z-values forming the surfaces of the DTM for the undisturbed plots occur significantly closer to the standardized base level (zero on the x-axis) than both the O/P washed and manually disturbed plots. However, if we assume that O/P washed sites were disturbed through HP/HW washing to a degree similar to our manual disturbance, this also suggested that these plots have recovered substantially since they were initially disturbed by the beach cleanup in 1989 and 1990. In Figure 12, we present the mean percentage of z-values extending varying distances above the standardized z = 0 (see Figure 2 ) elevation for a given sets of DTMs for all unwashed, O/P washed surveyed, and four plots manually disturbed during the study (KN5-1, KN5-2, KN208-1, and KN208-2). At the 50-mm slice, 31% of the z-values for the unwashed plots extend above the standardized zero elevation, compared to nearly 37% of the z-values for O/P washed plots. Finally, 73% of the z-values for the four disturbed plots extend above that elevation. This pattern continues for the remainder of the O/P washed and disturbed curves. This analysis demonstrates that a larger proportion of the z-values forming the surfaces of the DTM for the undisturbed plots occur significantly closer to the standardized base level (zero on the x-axis) than both the O/P washed and manually disturbed plots. However, if we assume that O/P washed sites were disturbed through HP/HW washing to a degree similar to our manual disturbance, this also suggested that these plots have recovered substantially since they were initially disturbed by the beach cleanup in 1989 and 1990.
section_name: 6.

section_name: Discussion

section_content: In its most simplistic form, cross-shore swash zone sediment transport is a function of the shear stress of wave swash and backwash over the bed and the critical shear stresses required to initiate movement of clasts. Once an armor has developed within coarse, surficial sediments, a stronger current is required to transport individual clasts, a phenomenon termed "structural strengthening" [12]. Organization of this armor layer further enhances the stability of the bed by maximizing clast-on-clast contacts, increasing the shear stress required to initiate movement of any one clast. Factors such as clast size, clast pivotability (relative projection of a clast above the bed, influenced by clast shape and orientation relative to surrounding clasts), degree of imbrication, and presence of a stable coarsegrained armor can alter this critical shear stress [5, 6]. The protrusion of clasts above the matrix bed level (increased pivotability) increases friction coefficients, thereby making the bed more erodible, whereas the most stable bed configuration is one fully armored, that is, one that has a high friction angle and low height above the bed where the flow is acting [39]. Thus, an armored bed that is organized is in a highly stable configuration. 

Our newly proposed metric for quantifying organization builds upon a wealth of work to quantify gravel bed structure or "roughness" in fluvial systems, including those
section_name: Utility of an Improved Organization Metric

section_content: In its most simplistic form, cross-shore swash zone sediment transport is a function of the shear stress of wave swash and backwash over the bed and the critical shear stresses required to initiate movement of clasts. Once an armor has developed within coarse, surficial sediments, a stronger current is required to transport individual clasts, a phenomenon termed "structural strengthening" [12]. Organization of this armor layer further enhances the stability of the bed by maximizing clast-on-clast contacts, increasing the shear stress required to initiate movement of any one clast. Factors such as clast size, clast pivotability (relative projection of a clast above the bed, influenced by clast shape and orientation relative to surrounding clasts), degree of imbrication, and presence of a stable coarse-grained armor can alter this critical shear stress [5, 6]. The protrusion of clasts above the matrix bed level (increased pivotability) increases friction coefficients, thereby making the bed more erodible, whereas the most stable bed configuration is one fully armored, that is, one that has a high friction angle and low height above the bed where the flow is acting [39]. Thus, an armored bed that is organized is in a highly stable configuration. 

Our newly proposed metric for quantifying organization builds upon a wealth of work to quantify gravel bed structure or "roughness" in fluvial systems, including those that apply similar digital mapping techniques [40] [41] [42]. However, in the wave-exposed (bi or multi-directional flow) gravel beach system, the consequences of exposure to water flow are quite different from what is reported in unidirectional flow streams or flumes. Notably, in the latter, exposure to downstream currents increases, rather than decreases, roughness [43]. Our approach is statistically robust, fast and easy to implement in field settings, and-because resulting OM values are independent of sediment grain size-appropriate for use across a wide assortment of gravel beaches characterized by clasts of widely varying sizes and shapes, and overall sediment compositions ranging from well-sorted to heterogeneous. As the case study shows, this approach can be employed on disparate intertidal heterogeneous gravel beaches that are subject to any manner of natural or anthropogenic disturbance. Further, OM serves as a proxy for surface clast stability in addition to clast orientation/angle. For example, a plot with partially buried clasts (and thus more muted plot microtopography; lower zbar and σ z values) would have a higher resulting value for OM than a plot with clasts of identical size and orientation, but in which those clasts sit entirely upon the surface of the bed matrix. This is likely one factor that accounts for the wide variation of OM at any specific value for abar (Figure 11 ). Likewise, OM values would be higher in sites with high degrees of imbrication, as the overlapping clasts themselves create the visible bed surface, minimizing microtopography above the z = 0 surface. The development of imbrication reflects extrinsic factors such as wave energy, as well as characteristics of the clasts themselves, such as clast shape: flatter clasts are more easily imbricated than clasts in which the c-axis approaches the dimension of the b-and a-axes. Because OM is normalized only by clast a-axes, it is insensitive to clast shape. However, imbrication is conspicuously absent from most PWS beaches [5] and was not observed in any of our study sites and cannot be addressed with these data. Thus, the response of OM to different degrees of imbrication remains an avenue for future exploration. 

The newly developed OM uses an off-the-shelf digital camera which, even combined with a proprietary software license and peripherals (estimated at <US$5000 total in 2021), cost far less than traditional ground-based LiDAR (3D laser scanner) systems. However, it would be difficult to adapt this system to subtidal or in-stream environments. Such applications would be particularly useful, for example, in assessing the effects of marine aggregate harvesting [44] or scallop dredging [45] on clast organization in heterogeneous gravel beds, which have been found to host rich sedentary megafaunal assemblages [45] [46] [47]. Because it is likely that clast organization is an important component in the recovery of disturbed megafaunal assemblages in coarse sediments [18], such knowledge would help in estimating the likely time required for recovery of the valuable targeted or untargeted living resources that are "harvested" during these operations. Recovery time is an important element in weighing the overall environmental cost/benefit of a proposed activity that is generally severely underestimated [18, 48, 49]. 

It is likely that newer technologies such as handheld (or now even phone-based) LiDAR and structure-from-motion could provide efficient alternatives to the PhotoModeler approach presented here. These have been demonstrated for use in similar-sized plots (e. g., [50] [51] [52], and at least the former could be applied in shallow-water environments. Further, expanding our approach with digital, in-situ grain-size analysis [53] would improve reliability of clast measurements used in OM for grain-size normalization (abar), and remove potential biases associated with Wolman Pebble Counts (e. g., [54]. However, care must be taken to ensure that clast-size measurements using such an approach are accurate and comprehensive. For example, in a case in which clasts are in a disorganized state (oriented on end, with 'a'-axes perpendicular to the beach slope), in-situ size measurements would be unable to capture the full length of the a-axes, as a substantial percentage of the length of the clast could be buried and not visible from the surface photogrammetrically, resulting in inaccurate values for abar.
section_name: Utility of an Improved Organization Metric

section_content: Our comparison of unwashed and O/P washed beaches in PWS reveals new insights into the long-term effects of HP/HW washing in both surficial clast organization and shallow subsurface stratigraphy. 

Overall, PWS beaches are coarser grained than heterogenous coarse beaches described elsewhere. For example, we found that gravel comprised 32% to nearly 100% (mean: ~75%) of surficial sediments, as compared with only 15-24% reported by Hayes et al. [2] for Prince William Sound beaches, Jennings and Shulmeister [3] in New Zealand, and Mason and Coates [4] from a literature review. At unwashed sites, we found that both Mz and percent gravel were higher in the upper and middle horizons and were lowest in the lower horizon (Table 5 ). Sand content was highest in the upper horizon. Silt/clay content was lowest in the middle level and highest in the lower horizon. None of these differences were significant. This is in contrast to observations of gravel stratification in beaches elsewhere, such as Wales [55] and Argentina [6]. A possible reason for this is that PWS beaches are geologically quite 'young', i. e., they were only established within the intertidal zone during the 1964 Great Alaska Earthquake, which elevated or lowered beaches considerably in western PWS [25, 56]. However, differences between unwashed and O/P washed beaches-in which the O/P washed beaches were significantly overall coarser and, at least within the upper and lower horizons, more gravel-dominated than unwashed beaches-importantly reveals that HP/HW beach washing altered the vertical structure of these PWS beaches down to at least 15 cm. Calculated values of OM were statistically independent of the broad range of grain sizes across our 37 sampled beaches in PWS, but the OM values were quite variable. This reflects the complexity of processes responsible for development of surficial organization in heterogeneous gravel sediments in the lower intertidal zone. For example, the responsiveness of clasts on a beach to reorientation into a more stable state (i. e., their pivotability) could be influenced by clast geologic composition (e. g., shape and specific density), differences in wave exposure or tidal currents), or the susceptibility of those clasts to reworking, as influenced by biological activity (e. g., cover of clasts by heavy algal cover, mussels, or barnacles, and excavation activities of predators). For example, the potential role of tidal forcing is suggested at IN32, which had the highest observed OM (13. 99 ). This site is located in a narrow channel protected from wave action, but exposed to strong, twice-daily, bidirectional tidal currents. Unfortunately, our small sample size of beaches unaffected by post-EVOS cleaning activities prohibits exploration of the role of these factors in sediment organization. This remains an opportunity for future studies and modeling efforts.
section_name: New Insights into Sedimentologic Structure of Heterogenous Gravel Beaches in Prince William Sound

section_content: Using western Prince William Sound (PWS), Alaska (USA) as a natural laboratory, we developed a new metric for quantifying the degree of organization of large clasts in surficial sediments in coarse, heterogenous beaches. This unique approach employs creation of three-dimensional digital terrain models from small study plots in the lower intertidal zone on a diverse series of heterogeneous gravel beaches in PWS. The resulting organization metric (OM) is independent of differences in clast size. Values varied widely among the beaches sampled, but OM is a useful tool for quantifying an important aspect of the organization on armored beaches. 

Overall, this study developed and field-tested a new approach for quantifying organization of heterogenous coarse sediments in diverse, terrestrial settings, including deserts. Although our particular approach relies on quantification of surface texture from photogrammetric analyses and field measurements of clast sizes, newer tools such as ground-based LiDAR and structure-from-motion could allow for faster application and more accurate measures of surface texture through digitized grain-size analysis. 

Finally, our case study using this new metric found that organization of the armored layer on beaches was significantly lower on those beaches which were observed or presumed to have been washed following the oil spill in PWS than on those which were left unwashed. Combined with observations of relatively coarser and more gravel-dominated stratigraphy to a depth of at least 15 cm in these washed beaches, this reveals that the physical structure and sedimentology of beaches that underwent post-spill washing remained significantly altered 21 years after the cleanup. This finding has important implications for beach management during and following major disturbance events (e. g., cleaning or dredging).
section_name: Conclusions

section_content: Acknowledgments: This is contribution 4088 of the Virginia Institute of Marine Science, William & Mary. The authors thank David and Annette Janka, operators of the R/V Auklet, and Terry Kennedy, Alaskan Wilderness Air, for valuable logistical assistance in data collection. Marcio Boechat Albernaz provided critical assistance in the field and conducted granulometric and statistical analyses for the bulk sediment samples in the laboratory. Emily Hein designed and ran MATLAB codes for analyzing photogrammetric data. We are also grateful to Matt Klymson and Dan Money of EOS ® Systems (Toronto, ON, Canada ) for assisting us in adapting the PhotoModeler Scanner system to our concept for measuring organization of heterogeneous sediment.

section_content: This research was funded by the NOAA/NMFS Exxon Valdez Oil Spill Office for the Exxon Valdez Oil Spill Trustee Council under Contract No. AB133F-10-CQ-0066 from Western Acquisition Division Contracts Branch to Littoral Ecological & Environmental Services. Institutional Review Board Statement: Not applicable.

section_content: The following supporting information can be downloaded at: https:// www. mdpi. com/article/10. 3390/jmse10040525/s1. Table S1 : Results of Granulometric Analysis.

section_name: Data Availability Statement:

section_content: All supporting data for this study are reported in the manuscript and/or supplemental materials. Raw data for sediments and photogrammetric analyses can be obtained from D. C. L. or C. J. H.
section_name: Data Availability Statement:

section_name: Conflicts of Interest:

section_content: The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
section_name: Conflicts of Interest: