Bathymetric Change Model of San Francisco Bay

Fregoso and others (2023) took bathymetric survey data collected from 1971 to 1990 (herein referred to as the 1980s bathymetric surface) and compared it to bathymetric survey data collected from 1999 to 2020 (herein referred to as the 2010s bathymetric surface). This comparison resulted in the bathymetric change surface (Fregoso and others, 2024) that allows for volume calculations of total sediment eroded from or deposited on the Bay floor during this period and serve as the foundation for our net sand volume change analyses. Our sand volume calculations are thus limited to within the spatial extent of the bathymetric change area shown in figure 2.

Sediment Cores

To quantify the volume of sand made available to the system (through bed erosion) and sediment sequestered (by deposition), we first determined how much sand is contained within the Bay floor and subsurface sediments. To do so, we created a 3D sand model to capture the variability in percent sand content both throughout the Bay and with depth beneath the Bay floor. The primary data for the creation of the sand model were sediment cores, predominantly gravity cores, but also box cores and vibracores, that were collected throughout the Bay in 1990 to 1991 (Anima and others, 2005), 2000 (Allison and others, 2003), 2006 (Yee and others, 2011), and 2016 (Takesue and others, 2021). Of the cores collected, 186 cores (ranging in length from 14 to 360 centimeters [cm]) were located within the 1980s to 2010s bathymetric change area in Fregoso and others (2024) which provided the volumes of sediment deposition and erosion used for this study (table 1, fig. 3). Few cores were collected in shallow, intertidal waters because of limitations in getting the boat into these areas and the necessity of having enough water for core collection to be successful.

Surface Sediment Samples

Examining the spatial distribution of the sediment cores revealed there were still large sections of the Bay without grain size information (fig. 3). To better constrain our 3D sand model, we decided to incorporate readily available bay-floor samples that have been analyzed for grain size. Bay floor samples are not ideal for creating a 3D sand model because subsurface grain size may differ from grain size at the surface of the bay floor. In fact, the cores used in this study proved that there are definite variations in sediment size within the length of a core, but given the lack of available core data, and following suit with a study by Alonso and others (2021), we applied the assumption that subsurface grain size is the same as the surface size in locations where that was the only available data.

Two readily available sources of quality floor-surface grain-size data were used in this study. An extensive sediment provenance study by Barnard and others (2013) examined surface sediments (top 10 cm) from a total of 339 grab samples collected in 2010 to 2012, 209 of which were within our study area. These samples reported grain size as percentages by size class, calculated by weight. All 209 of the samples from the sediment provenance study contained some percentage of sand. Using these samples alone would have biased the models towards sand, as mud-only areas were not adequately defined to aid in the interpolation between cores and shallow waters near the shoreline.

For creating a sand budget, modeling areas with an absence of sand are just as important as areas with sand. The Bay Regional Monitoring Program (RMP) for Water Quality database of long-term sampling efforts in the Bay (San Francisco Estuary Institute, 2023) was used to obtain samples denoting areas without sand. Through the RMP website (https://www.sfei.org/programs/bay-regional-monitoring-program) it is possible to download sample data, by analyte and (or) fraction size, for sediment samples collected during a specified time-period. Sampling for this study ranged from 1993 to 2018, with many of the same locations repeatedly sampled throughout the years. The RMP sample database does not give a definitive weight percent of each grain size class (for example, sand, mud, and so forth) for every sample, so it cannot be used like the Barnard and others (2013) samples to constrain the exact sand content without more extensive, time-consuming work combing through the dataset to organize the data. However, the RMP database does provide a descriptive designator like the ones that Folk (1974) used (for example, fine silty mud, muddy sand, and so forth) that allows a range of sand contents for the same descriptor for every sample. The initial download of RMP sediment sample data, of size classes for sediment less than 2.0 mm (sand, clay, silt, fines [silt+clay]), returned more than 7,000 samples. Of those, 498 samples were ultimately used in this study (refer to the “Methods” section). Figure 4 shows the distribution of the surface samples that were used in relation to the core locations.

Areas of Human Disturbance

Within each of San Francisco Bay’s subembayments, there are regular activities that alter or maintain the depth of the Bay floor. These projects include dredging of channels for navigation, disposal sites for dredged and other sources of sediment, and shell and sand mining. Figure 5 shows areas of such activities within the Bay. These areas are delineated by a combination of information on current activities: U.S. Army Corps of Engineers (USACE, 2023) hydrographic surveys of the Federal navigation channels and dredge material placement sites, current sand mine leases (California State Lands Commission, 2012), and unnatural patterns of bathymetric change from the 1980s to 2010s (Fregoso and others, 2023). These areas, due to their heavily managed states, are not indicative of natural influences. Because of this, any grain size data within these areas have been excluded from sand content interpolation with the only anthropogenic exception being the sand mining areas. These areas are being mined because they contain naturally occurring sand beds. Therefore, the presence of sand in these areas are included in the interpolation process through the addition of simulated cores within the sand mine areas with a sand percent value of 100 assigned for the entire length of the core.

Areas of Oyster Shell

South San Francisco Bay is a historical source of native Ostrea lurida (Carpenter, 1864) (Olympia oysters) that have been collected, first by the local indigenous people, and then mined for more modern uses beginning as early as the 1890s, with more industrial efforts starting in the 1920s and still occurring today (Hart, 1966). Oyster and shell hash have a known presence in many areas of the Bay and can be inadvertently classified as sand through simple grain size measurements that do not distinguish between organic and non-organic particles. To prevent conflating sand-sized shell hash with mineral sand grains, the large shell beds in south Bay, identified with an acoustic seabed classification survey in 2005 (acoustic class 10 in sub-appendix B in Foxgrover and others, 2007) were excluded from sand content interpolation (fig. 5) and volume estimates provided as separate results.

Interpreting and Preparing Core Data

The production of the 3D sand models was a multipart effort that began with the analysis of core data in the form of core logs. Core logs are detailed handwritten descriptions of cores recorded by researchers noting physical changes in the cores including measured changes and thicknesses of the sediment layers within the core, organics such as roots and wood, presence and type of shells, and color changes. The core logs are hand drawn representations of the core using standardized geologic symbols with written descriptions detailing, for example, more complicated areas, marking the specific depths of changes in the core, and the start of a new layer of sediment denoting a distinct sedimentary unit within the core. These logs were provided in the form of scanned paper logs, or in cases of published cores as graphicly illustrated logs which have been produced for a report. Some of the cores also had X-radiographs done, and measurements from sediment grain-size analyses at specified locations within the core. Examples of core logs are provided in figures 1.1 and 1.2 of appendix 1.

To determine the percent sand content throughout the length of each core, the available data were examined, including a small number of grain size analyses and graphical logs contained in reports (Woodrow and others, 2010, 2014, 2017a, b). All cores were measured in centimeters. Where sand was described in the written portion of the logs and (or) drawn in the visual depiction of the core, percent sand was estimated based upon the information provided. The illustration of the core helped contextualize the distribution and thickness of sand in the intervals defined by the core log. Several logs included a written percent sand estimate for intervals, which were then directly transferred to a spreadsheet for this effort (Fregoso and others, 2025). Where there was no percent sand estimate, various methods were used to estimate sand presence in the interval. If an interval contained sand laminae, then they were counted and estimated for thickness from the visual depiction of the core and associated descriptions to determine the amount of sand in the interval. If an interval contained sand interbeds or lenticular bedding, a similar process was used depending on the thickness of each bed. Where there were sand-filled burrows, an estimation was made based on the frequency and distribution of the burrows and their apparent thickness drawn on the log. If an interval was described as “muddy sand,” “sandy mud,” or any other description in which sand was either a matrix or a secondary grain size in any interval and not clearly distinguished visually on the log (that is, sand was noted in the description but not visually marked in any identifiable way on the drawn core), a best estimation was made based on the description. Percent sand content for core intervals was calculated by dividing the total amount of sand in centimeters from the length of the interval in centimeters and multiplying by 100.

The core data analysis was divided into separate Excel files for each subembayment. Each core had its own sheet that reported these sand percent estimations as well as minimums and maximums for each estimation. When estimates for sand were based on visual representations and were given less precise numbers, such as 40 percent, minimum and maximums were assigned as plus or minus 10 percent. When estimates were determined by counting and measuring specific intervals and to calculate a specific percent, such as 5 cm out of 20 cm for 25 percent sand, the minimum and max values were plus or minus 5 percent. The reports also include general observations of sand. For example, details are given on grain size, laminae, bedding, burrows, and difficulty of making such estimations depending on available detail on the paper logs (table 2).

To create the 3D sand models, data from the core data spreadsheets were reformatted for use in geographic information systems (GIS). The sand model depicts percent sand content in 10-cm-thick sections throughout the length of the cores. To do this, the data shown in table 2 were adjusted to be displayed in regularly spaced, 10-cm intervals referenced to the top of the core. If the top of the cores had gaps between the top of the core liner and the beginning of the actual sediment core, the distance of the gap was subtracted from the raw depth intervals. In the case of the core in table 2, for the first section interval, 2 cm was moved to zero (0) cm, and 41 cm became 39 cm. All intervals had to be adjusted by 2 cm as shown in table 3.

Once the core was vertically adjusted, a weighted average was used to determine the sand percentage for each 10-cm interval of the core. For the section 30–40, 30–39 has a value of 5.13 percent, where 39–40 has value of 0 percent. So, the value for 30–40 is

There were some instances where cores had two different logs from examinations in different years, which meant two separate core analyses. When there was a conflict between the two analyses, the midpoint was taken, and the “min” of the minimums, and “max” of the maximums were used. If one log showed sand, and the other did not, the one with sand was used. Where there were grain-size measurements, which never represented more than 3 cm in an interval, those values and locations within an interval were used in the weighted average calculation. Grain-size measurements were considered absolute, with no difference in min or max values. Although percent sand amount in grain-size measurements was calculated by weight and the estimates from core logs is percent sand by volume, this affected less than 1 percent of the data, the 0.86 meter (m; 3 feet [ft]) that were analyzed for grain size in a total of more than 265 m (869 ft) of core length that had sand present. The lack of a conversion from sand percent by weight to sand percent by volume does not create a significant error compared to the other uncertainties in the analysis (refer to the “Uncertainty in Sand Volumes” section). The sand content values for all 186 cores were compiled into a single spreadsheet (Fregoso and others, 2025); each core has a row in the table with its associated x-y location data and a column containing sand percentages for each 10-cm interval.

Core collections occurred in 1990 to 1991, 2000, 2006, and 2016. The core collection dates were not coincident with the bathymetric survey dates that constrain the volume change calculations. To accommodate this, the cores were vertically shifted to account for changes in bed levels that occurred between the time the cores were collected, and when the elevation of the Bay floor was surveyed in the 1980s by using the yearly rate of change in the seafloor from the 1980s to the 2010s. To account for these differences, the vertical reference of the top of a core was shifted to coincide with an interpolated floor-surface for the year of its collection. This was calculated using the following variables:where

For core 90–134 (located in central Bay, Anima and others, 2005), b=1979, s=41, c=1990, and chg=1.07 m.

This example is for a core where there was deposition that occurred after the core was collected. The equation results in 0.3 m, which indicates that the entire core needs to be shifted upward 30 cm to level the core to the 1979 bathymetric surface reference plane for central Bay (fig. 1). This means core 90–134 with an original length of 238 cm now has 30 cm above the 1980s baseline, and 210 cm below (table 5). Had the above equation calculated negative 0.3, indicating erosion of sediment occurring after the date the core was collected, all the values for the entire core would have been shifted downwards 30 cm. The values from 0 to 30 cm would be assigned the same value as the top of the core. In this case, the values for the first row in table 4, would be repeated up to and including the 30–40 cm interval; the bottom of the core would now end at 270 cm below the surface.

The final adjustment for the core data prior to generating the 3D sand model was to convert the interval data located from separate worksheets to one table that could be easily imported to ArcGIS Pro 3.3 (Esri, Redlands, California) for geospatial analyses. This required the information for each interval to be converted from rows to columns for each core. Each core has a row in the table with its associated x-y location data, and a column for each 10-cm interval containing percent sand content (Fregoso and others 2024, 2025). These new tables can then be used to plot the core locations on a map with all its core interval data (table 6) and generate raster surfaces of percent sand content for each 10-cm-thick horizon.

Preparing Surface Sediment Samples

Including the Bay floor sediment samples in the grid interpolations required checking for conflicts at locations with the cores, filtering the RMP samples to select mud only samples, and converting sand percent by weight to sand percent by volume to match the core data. To minimize conflicts between two different data types while ensuring the priority of the sediment core data during the grid interpolation process, all surface samples located within 300 m of a core were removed.

The RMP database of 7,000 samples was analyzed to identify mud-only cores. First, all samples were plotted on a map of the study area, and samples collected outside of the 1980s to 2010s bathymetric change area were removed. Next, samples that were within 100 m of a mud-only sample were selected and examined to see if any samples contained sand. If a location revealed a sample containing any sand, all samples within that 100-m radius were removed, including mud-only samples. Finally, any remaining samples with sand content were removed. This selection criteria resulted in 498 data points with a high confidence of representing areas of mud, or zero (0) percent sand.

The sand samples from the sediment providence study (Barnard and others, 2013) reported percent sand content calculated by weight. For these samples to be directly comparable with the percentages from the core analyses, these percentages needed to be converted to percent by volume. Following Alonso and others (2021, their equations 3, 5–7), the sand percentages in the surface samples were converted from a weight-derived metric to a volumetric one using a dry bulk density value of 1,530 kilograms per cubic meter for 100 percent sand (Mulder, 1995, from Alonso and others, 2021).

Creating a 3D Sand Bed Model

To calculate the amount of sand stored or released into the system, including potential uncertainty in results stemming from variability in how the sand content layers were modeled, a total of four 3D sand models were generated for each subembayment, each with 100×100-m cells representing a thickness of 10 cm. First, an all-inclusive 3D model was generated from the best estimate of sand content from cores, sand mine simulated cores, and surface samples. Then, for use in assessing variability in uncertainty, three sand models were generated from only the cores and sand mine simulated cores, varying only the range in values reported in the core logs (1) using values for best estimate, (2) using minimum, and (3) maximum sand percentage values.

To convert between point measurements and grids of percent volume of sand content, different methods of interpolation were explored: Inverse Distance Weighted (IDW), natural neighbor, spline, and topo to raster. For this study, the tool that easily allowed multiple input files for interpolation without the need to combine all data into one file was selected, the topo to raster tool (Esri, 2023). The topo to raster tool allows us to add or subtract inputs, cores and later surface sample data, generate different versions, adjust the controls on interpolation, input minimum and maximum values, set allowable tolerances, and vary the degree of smoothing to create a model that best fit the data.

The 3D sand model’s size is strictly defined as the boundaries, both in the horizontal and vertical, of the bathymetric change grids (Fregoso and others, 2024). For instance, in the central Bay sand model, the maximum deposition is 1,320 cm and the maximum erosion is 1,660 cm. The vertically referenced cores range from 60 cm above the 1980s reference plane for deposition, to 240 cm below the reference plane for erosion. Because the cores do not extend into the full extent of the depositional or erosional layer, the sand percentages at the top and bottom of the cores are extended to the maximum bathymetric change values. By extending the values, this allows for all core locations to contribute data for interpolation purposes for each layer in the 3D sand model. Layers generated for the minimum and maximum range of the cores within each subembayment are applied to all sediment volumes either shallower than or deeper than the extent of the core data.

To help assess variability due to the relatively sparse core data and its effect upon sand volume estimates, we carried out additional analysis limiting calculations to areas of the Bay well constrained by sediment core data. To delineate such areas, we compared the surface layers of the best estimate of the cores-only models with the surface layer of the surface sediment sample layers. We did this by taking the top 10 cm layer from the adjusted and (or) leveled cores from the best estimate model (fig. 6) and compared it to the top 10 cm layer of the all-inclusive model (fig. 7). By using just this surface, the comparison was between the actual data produced from the core, no extrapolation, to data produced from the cores combined with the surface samples. The two models were differenced, subtracting one from the other for each embayment. This differenced surface reveals where the core-only interpolations were similar to the interpolations supported by greater data density, core and surface sample data. Areas of differences less than 10 percent (fig. 8) define an area of higher confidence, where the interpolation in the core-only layers of the 3D sand model are best constrained by actual core data. Calculations limited to these high confidence areas (fig. 8) are presented within our uncertainty analyses (refer to the “Uncertainty” section).

Calculating Sand Volumes

The technique for calculating the volume of sand in areas of sediment deposition or erosion during the 1980s to 2010s is based on the methodology used by Foxgrover and others (2019) to estimate legacy mercury remobilization in Alviso Slough in the lower south San Francisco Bay. For each subembayment, the bathymetric change grids from Fregoso and others (2024) were parsed out into separate grids for areas of sediment deposition versus areas of erosion to calculate the volume of sand stored (deposited) or released (eroded) during the 1980s to 2010s. The deposition and erosion grids were then sliced into 10-cm-thick horizons, and the volume of sediment calculated on a cell-by-cell basis by multiplying the horizon thickness by the raster cell size (25 m×25 m for Suisun, San Pablo, and central San Francisco Bays, and 50 m×50 m for south Bay). The percent sand content grids from the 3D sand model were resampled to match the cell size of the deposition and erosion grids, 25 m or 50 m. Then, for each 10-cm-thick horizon, the sediment volume layers were multiplied by the corresponding percent sand content grids from the 3D sand model to calculate volume of sand. The individual cells were summed for each horizon and the total sediment volumes and sand volumes exported to a text file. Values were calculated for the all-inclusive 3D model, and the 10 percent clipped versions of the best estimate, minimum and maximum models.

For areas designated as affected by human disturbance, including sand mining, dredging and disposal areas (fig. 5) sand volumes were calculated using average sand content values, provided by subembayment in McKee and others (2023) and bathymetric change from Fregoso and others (2024). Sand mining areas were assumed to contain 100 percent sand, disposal in central Bay was assigned 16 percent, and disposal in San Pablo Bay assigned 45 percent. Sand percent for dredging areas differed by subembayment: Suisun Bay and Carquinez Strait were assigned 7 percent, San Pablo Bay 58 percent, central Bay 30 percent, and upper south Bay 7 percent. The oyster shell hash areas in south Bay used the same methods as above to determine sand content, but with an assumed sand content of 10 percent.

Sand Volume Change in the Bed Excluding Areas Affected by Human Activities and South Bay Oyster Beds

The 767 km² of study area that was not obviously affected by anthropogenic activities underwent a greater amount of deposition than erosion of bed sediments. However, there was a net loss in the amount of sand stored in the bed equating to −2.64 Mm³ (table 7). This area excluded mining, dredging, disposal areas, and the south Bay area containing oyster shell deposits. A negative change in net sand volume in bed sediments indicates that the bed is a source of sand—more sand was eroded from the bed than deposited. The largest volume of sand, approximately 4.36 Mm³, lost from the bed through erosion in the study area was in San Pablo Bay. Erosion in Suisun Bay, Carquinez Strait, and upper south Bay resulted in a net decrease in the amount of sand stored in the bed by about 5.42 Mm³. In contrast, about 6.40 Mm³ of sand was added to the bed in central Bay through deposition, in the areas not mined for sediment. In lower south Bay outside the area of the oyster shell deposits, 0.74 of Mm³ of sand was added to the bed by deposition (table 7).

Both the total sediment volumes and the sand content of sediment eroded and deposited affected the net change in sand storage in bed sediments. For Suisun Bay and Carquinez Strait, both the greater volume of erosion (23.55 Mm³) than deposition (8.96 Mm³) and the higher sand content in eroded sediment (91.8 percent) than in deposited sediment (44.2 percent) resulted in a net loss of sand from the bed (table 7). In contrast, for San Pablo Bay, the net loss of sand from the bed was driven by the greater volume of erosion (41.97 Mm³) than deposition (32.2 Mm³) and not the sand content.

The rates of net change in sand volume (in millions of cubic meters per year [Mm³/yr]) also were calculated to account for the variation in subembayment-averaged time spans (20–35 years) between the 1980s and 2010s bathymetric survey collection dates. San Pablo Bay had the highest net rate of sand erosion (−0.14 Mm³/yr) for a span of 31 years. Upper south Bay (29–year span) and Suisun Bay and Carquinez Strait (20-year span) both had erosion rates of −0.11 Mm³/yr of sand from their beds. Central Bay, with the longest time span between surveys of 35 years, had the highest rate of deposition, 0.18 Mm³/yr, and lower south Bay with a span of 28 years, rate of sand deposition was, 0.03 Mm³/yr.

Sand Volume Change in the Bed in Areas Directly Affected by Human Activities

The net sand volume loss from the bed in areas directly affected by human activities (mining, dredging, disposal) (fig. 5), encompassing 45 km², was about −14.06 Mm³ (table 8). Both the volume of sand loss and the rate at which it was lost, −0.46 Mm³/yr, are about twice that of the area not directly affected by human activities, while only being 13 percent of the total study area. About three-quarters of the loss was in sand mining areas, with most (about 70 percent) occurring in Central Bay. The sand volume loss in sand mining areas reported here is based upon how much the elevation of the bed lowered between two discrete points in time. Although related to the volume of sand removed from the bed by sand mining, which varied during our study period and was estimated by Perry and others (2015) to average −1.1 Mm³/yr removed from 2005 to 2014, bathymetric change within mining areas also includes sediment deposition and erosion from natural processes. From the 1980s to 2010s, the average lowering of the bed in the mining areas of central and Suisun Bays were −1.0 m and −0.8 m, with rates of sand loss from the bed of −0.25 Mm³/yr and −0.11 Mm³/yr, respectively (table 8). Roelvink and Elias (2023), studied bathymetric change in San Francisco Bay from 2008 to 2019 and found a net decrease in sand storage in the bed from bed lowering in sand mining areas in both central and Suisun Bays, except for one lease area in central Bay during specific time periods.

Net sand loss from the bed in dredging areas was about 4.1 Mm³ and occurred at rates ranging from less than −0.01 to −0.06 Mm³/yr (table 8). Net sand gain to the bed in disposal areas is small, about 1.06 Mm³ in total, compared to a net loss of about 15.1 Mm³ from sand mining and dredging areas.

Sand Volume Change in the Bed in South Bay Oyster Shell Deposit Areas

The bed of much of the upper south Bay (65.3 km²) and parts of lower south Bay (0.3 km²) are composed of a mixture of mud and native oyster shells (Foxgrover and others, 2004). Visual observation reveals that the sand-sized particles in these oyster shell deposits are overwhelmingly broken shells, but mineral grains that are sand sized are also present. If grain size analyses were carried out on samples after dissolving the shell material with acid, the mineral grain sand content would be measured. Unfortunately, analytical steps did not include acid dissolution to remove shell material, so the resulting sand content values are for a mixture of shell and mineral grains. If we assume that 10 percent of the sediment eroded and deposited is sand-sized mineral grains, 0.55 Mm³ of sand (defined as mineral and not biogenic material) was eroded from the bed from the 1980s to 2010s at a rate of 0.02 Mm³/yr (table 9). The contribution of oyster shell deposit areas to the total loss of sand from the bed in the Bay is small compared to other areas and because the assumption of 10 percent sand content is not supported by data, we do not include the hypothetical 0.55 Mm³ sand loss in our total estimate.

Uncertainty in Sand Volume Change

Uncertainty in our estimates of the change in sand storage in the bed of San Francisco Bay derives from both inexact knowledge of sand content in the bed and uncertainty in bathymetric surveys, which affects the volume estimates of sediment deposition and erosion between surveys. The uncertainty in estimates in sand content has three primary sources: (1) interpretation of core logs, (2) differences in sand content in the upper 10 cm of cores and surface samples collected in the same areas at different times, and (3) the assumption of constant sand content with depth below the surface at surface sample locations. Analysis of the effect of uncertainty in bathymetric surveys used in this study on sediment volume change is detailed in Fregoso and others (2023). Here we present a modification of that analysis for estimating uncertainty in sand volume change.

Uncertainty from Interpretation of Core Logs

Estimating the sand content in cores from graphical logs is fundamentally subjective. To help quantify the effect of this subjectivity, we developed sand models using the full range of sand content values provided by the researchers that initially analyzed the cores and created the logs. Table 10 contains volumes of sand content derived from reported best estimates of sand content, as well as estimates of the minimum and maximum sand content of core-only 3D models. These estimates were made for the portions of the Bay within the defined 10 percent high confidence clip (fig. 8) where the density of cores is high enough to constrain spatial interpolations between core locations without the additional use of surface sediment samples. In total, the high confidence area covers about 527 km² or about 69 percent of the study area and excludes the south Bay oyster shell beds and areas directly affected by human activities.

The effect of sand content estimates on the volumes of sand deposition, erosion, and the change in net sand volume in the bed for areas with a high density of cores in San Francisco Bay is presented in table 10. The range of net sand volumes stored in bed sediments using the minimum and maximum estimates of sand content in cores, made for each depth interval identified in the core, is −4.17 to −2.68 Mm³. The minimum estimate resulted in a net sand volume loss of 0.92 Mm³ less than the best estimate. The maximum estimate resulted in a net sand volume loss of 0.57 Mm³ more than the best estimate.

Uncertainty from Augmenting Core Data with Surface Sediment Sample Data

Because the density of cores was less than optimal, we augmented core data with floor-surface sample data. There are two potential sources of uncertainty associated with augmenting core data with surface sediment sample data: (1) differing sand content in surface sediment samples and the upper 10-cm layer of cores, and (2) variation of sand content with depth below the Bay floor. Both sources of uncertainty are difficult to assess.

We evaluated the effect of using surface sample data on net sand volume change in the bed by creating 3D models with (fig. 7) and without (fig. 6) augmenting the core sand contents with surface-sample sand contents in the less than 10 percent difference areas as defined in the “Methods” section. The 3D model including surface samples followed Alonso and others (2021) approach of extrapolating the surface sand content to all intervals below the surface—a constant sand content with depth.

The effect of augmenting core data with surface sample data is shown in table 11. The area of comparison is the same 527 km² as for the evaluation of uncertainty associated with interpretation of sand content from core logs. Adding surface sample sand contents decreased the net loss of sand from the bed by 0.39 Mm³ from 3.60 to 3.21 Mm³, which is an 11 percent reduction. The reduction in loss was caused by a higher sand content of surface samples in depositional areas than in erosional areas resulting in a greater increase in average sand content in depositional areas (1.2 percent increase vs. 1.1 percent increase in erosional areas).

The second uncertainty in using surface samples in addition to cores in the 3D model of sand content is that we assumed uniform sand content with depth following Alonso and others (2021). This assumption is incorrect for most individual cores used in this study but may be acceptable on average; sand content variations with depth below the bay floor may be sometimes less and sometimes more than at the surface. The average difference between the sand content for all 10-cm intervals below the 0–10 cm top interval in the 186 cores used in this study was 1.5 percent, indicating that, on average, the sediment below the surface interval is sandier. Table 12 gives the average sand content as about 29.8 percent for the areas of the 3D model with core data. Therefore, a 1.5 percent increase in sand content below the upper 10 cm of the cores equates to a 5 percent relative increase. The assumption of uniform sand content, on average, does not result in a large uncertainty.

Uncertainty Associated with Bathymetric Change

Fregoso and others (2023) analyzed volume uncertainty from bathymetric change, which is for all sediment sizes (sand, mud, and so forth). The effect of uncertainty in bathymetric change on our estimates of net sand volume change in the bed of San Francisco Bay is similar to, but not the same as, uncertainty in volume change for all sediment. The difference is that the volume of net sand change is a fraction of the volume of net total sediment change. This fraction is the ratio of the sand content of the bed sediment. For example, if the sand content of the bed is 30 percent, the uncertainty in the volume of net sand change is 30 percent of the uncertainty in volume of net total sediment change. Refer to Fregoso and others (2023) for uncertainty in bathymetric change. Here, we present key concepts and extend the analysis to uncertainty in sand volume change in the bed.

Two types of uncertainty are associated with bathymetric change: random error and systematic biases (USACE, 2013; Anderson, 2019). Random error is associated with sounding inaccuracy or noise and is generally randomly distributed in space (Adams, 1942; Shalowitz, 1964; Sallenger and others, 1975). Sounding errors are primarily low magnitude, lie both above and below the true value, and with enough data points cancel to negligible levels when averaged or summed as is the case in calculation of volume change (Anderson, 2019). Systematic biases, however, can enter surveys through various sources, including differences in horizontal or vertical datums (reference points), measurement inconsistencies, or during the creation of bathymetric grids. In comparing bathymetry from two periods for an area with many soundings, random errors cancel out and therefore do not significantly affect our estimates of deposition or erosion. However, systematic biases introduce false offsets between surveys and must be accounted for in bathymetric change estimates.

To minimize systematic error, Fregoso and others (2023) made efforts to understand and account for any differences in horizontal and vertical datums to not influence bathymetric change estimates. However, other sources of systematic errors are still present in our analysis, some of which are known and can be removed and others which are either unknown or difficult to remove; refer to Byrnes and others (2002) for a discussion of sources of systematic survey errors.

Esposito (2016) carried out an error assessment for the U.S. Department of Energy’s Office of Coastal Management (OCM) surveys, which are the majority of the 2010s surveys used to calculate the net change in sand volume in the bed of the Bay. Esposito compared independent depth values collected by different sonars and survey vessels in areas of overlap between adjacent survey blocks. The mean difference of the 1-m grid cells in areas of overlap was 0.019 m and the standard deviation, a representation of random error, was 0.144 m (Esposito, 2016). The mean difference in overlapping grid cells was interpreted by Fregoso and others (2023) to be caused by a systematic bias in one or both surveys mapping the overlap area.

Because the assessment of systematic error is not perfect and does not include all possible sources, and because systematic error likely varied for different parts of the Bay, a single value that we are confident characterizes the uncertainty from systematic error does not exist. However, based on the 0.019 m error found by Esposito (2016) in his error assessment of the OCM surfaces and our own similar finding of 2 cm or less of vertical uncertainty in USGS bathymetric surveys (Foxgrover and others, 2011), Fregoso and others (2023) found it reasonable that the minimum uncertainty of this change surface would be ±4 cm, the sum of 2 cm of uncertainty for the 2010s surveys and an estimated “best case” 2 cm uncertainty for the 1980s surveys. They presented volume change associated with two different possible levels of uncertainty: a minimal value of 4 cm, and a more conservative value of 8 cm that assumes uncertainty in the 1980s surveys is 6 cm. Depending on the application, they left consideration of uncertainty level to the user. Fregoso and others (2023) specify that systematic biases of multiple surveys could be random and could cancel out just as random sounding errors do when there are enough points. If this occurred, then uncertainty would be less than ±4 cm.

Table 12 shows the net sand volume change associated with 4- and 8-cm uncertainties in bathymetric change. The uncertainty in change in net sand volume in the bed associated with uncertainty in bathymetric change is calculated as either 0.04 or 0.08 m times the area of interest (in square kilometers) multiplied by the average sand content for that area. Table 12 reports uncertainties for areas without obvious human activities (table 7), areas with human activities (table 8), and oyster shell areas (table 9). Overall, for 4- and 8-cm uncertainties in bathymetric change, the uncertainties in net sand volume change in the bed in the entire study area are about 9.1 and 18.2 Mm³, respectively. The large area results in a large uncertainty. For reference, the Bay, excluding areas affected by human activities and oyster shell deposits, had a net loss of about 2.64 Mm³ of sand from the bed from the 1980s to 2010s (table 7), which is a smaller volume of sand than the 9.1 Mm³ associated with 4 cm of uncertainty in bathymetric change. The areas affected by human activities had a net loss of about 14.06 Mm³ (table 8) of sand from the bed from the 1980s to 2010s, which is more than 20 times the approximate 0.6 Mm³ associated with uncertainty in bathymetric change at the 4-cm level. The small area results in a small uncertainty.

Summary of Uncertainties in Sand Volume Change

We explored four sources of uncertainty in estimates of the net change in sand volume in the bed of San Francisco Bay: (1) uncertainty associated with estimates in sand contents of cores, (2) uncertainty from augmenting core sand content data with surface sample sand content data, (3) the assumption of uniform sand content with depth we made to use surface samples in a 3D model, and (4) the uncertainty associated with potential systematic errors in the bathymetric change analysis.

The three uncertainties associated with how sand contents were treated in the 3D model are all relatively small compared to the net change in sand volume in the Bay bed, which lost about 17.3 Mm³ combined (including areas of human activities and oyster shell beds) from the 1980s to 2010s. The uncertainty in net volume change associated with interpretation of sand contents from core logs, half of the range between minimum and maximum sand content estimate (table 10), is about 0.8 Mm³ for the areas with cores, which if scaled up to the entire study area is about 1.4 Mm³. The uncertainty in volume change associated with the addition of surface-sample sand contents (table 11) is about 0.4 Mm³ for the areas with cores, which if scaled up to the entire study area is about 0.7 Mm³.

Calculations were not made to determine the effects of potentially underestimating or overestimating sand content. However, on average, if sand content is 1.5 percent higher below the upper 10 cm of the core for a change in net sand volume in the bed, a comparison to the average sand content of about 30 percent indicates that the assumption of uniform sand content with depth would not have a large effect on net sand volume change in the bed. Together, the treatment of sand content introduced an uncertainty of about 3 Mm³ into the estimate of net change in sand volume in the bed.

In contrast, potential uncertainty associated with a systematic error of 4 cm in bathymetric surveys may result in 10.0 Mm³ of uncertainty in estimates of the net change in sand volume in the bed of the study area (table 12). This is about 13 times the uncertainty associated with sand content uncertainties (table 10). There is a 10.0 Mm³ uncertainty for 4 cm of bathymetric change uncertainty and if the uncertainty in bathymetric change is 8 cm the uncertainty in net change in sand volume in the bed doubles to 20.0 Mm³. However, if systematic errors in the surveys are not correlated, then the uncertainty could be less than 10 Mm³.

Limitations of This Study

This is the first study to our knowledge to attempt to estimate the change in net sand volume in the bed of San Francisco Bay. During our research, we discovered limitations in the study, including a lack of sufficient data. First, there were no bathymetric change data for parts of the Bay, most notably within Suisun Bay, which did not allow calculation of an entire Bay net-change in sand in the bed. In addition, surveys were conducted at different times for different parts of the Bay, so the picture is not synoptic for a changing system, which the Bay is. Although we assume that rates of bathymetry change are for the same period across the entire Bay, in reality they are not. Another data limitation was the paucity of cores in the Bay. In total, there were 186 cores available in the Bay, which is an average density of about one core every 5 km². For comparison, Dam and others (2022) based their study of the Western Scheldt in the Netherlands on 23,000 boreholes, which is an average density of about one core for every 0.7 km² that supports a 3D model with 100 m×100 m grid cells. The Alonso and others (2021) study of the western Dutch Wadden Sea used 4,800 surface samples to characterize sand content for an area smaller than the Bay. In this study, 209 surface samples from a sediment provenance study (Barnard and others, 2013) had grain size data detailed enough for our use, with another 498 RMP samples (San Francisco Estuary Institute, 2023) that we used to constrain the 3D sand model in areas where sand is absent.

The biggest limitation, however, was not data. The inability to unequivocally quantify the amount of uncertainty associated with the bathymetric change component of the change in net sand volume in the bed of San Francisco Bay calculation, limits interpretation of the results. For areas not affected by human activities, if a conservative estimate of bathymetric uncertainty is used, the uncertainty is similar to our estimates of net change in sand volume in the bed of the Bay. The conservative estimate is a worst-case scenario; the likely systematic error is probably much smaller (Fregoso and others, 2023).

Future Considerations

There are several next steps in this research that could make it more useful in developing sand budgets and improving the understanding of how sand is eroded, deposited, and moves in San Francisco Bay. A first possible step is to carry out bathymetric surveys in areas that were not surveyed in the 2010s, especially Suisun Bay. Analysis of such a survey would allow a more complete picture of the change in sand storage in the bed of the Bay. An analysis of gradients in sand content in the bed could inform where additional cores are needed to improve the 3D model of sand content in the bed developed for this study. An analysis exploring the effects of removing cores used in development of the 3D sand model could shed light on errors introduced by interpolation between cores.

The data developed in this study can also be used to calibrate and validate numerical hydrodynamic/sediment transport/geomorphic change models to both improve a sand budget and to better understand sand transport in the Bay. Such models can explore how sea level rise, climate change, and change in human activities could alter the sand system. Recent research that uses coupled hydrodynamic/sediment transport/geomorphic change numerical models to forecast the effects of sea level rise on geomorphic change has already started (for example, Elmilady and others, 2019; van der Wegen and others, 2019).

The application of the updated version of the SedTrails tool (Elias and Pearson, 2020) to the data developed in this study is a potential future step. SedTrails post-processes numerical sediment transport model outputs and can determine sand transport pathways and connections between different parts of a San Francisco estuary system by tracking the paths and number of sand particles travelling between a start and end point and how long the sand particles remain at the end point. To date, SedTrails has only been applied to two-dimensional model outputs. In the future, it could be possible to apply SedTrails to outputs from a 3D numerical model to account for 3D circulation (and sediment transport) that are present in parts of the Bay.

As mentioned in the “Limitations of This Study” section, an important future addition to this research would be to improve the estimates of uncertainty. Although not an easy task, more precise uncertainty estimates will increase the utility of this research. A component of this research would be to determine if the uncertainty is aleatory or epistemic, both of which are present in our study. Aleatory uncertainty, which is also referred to as stochastic uncertainty, variability, or irreducible uncertainty, derives from the natural randomness in processes and is typically characterized by probability distributions. Epistemic uncertainty, which Roy and Oberkampf (2011) refer to as reducible uncertainty or ignorance uncertainty, can be reduced or theoretically eliminated by increased understanding of the system. Research increasing the understanding of the aleatory and epistemic contributions to bathymetric uncertainty for surveys of the Bay, the largest source of uncertainty in our study, could improve the usefulness of the results. Because aleatory uncertainty in bathymetric surveys—random error in soundings—tends to cancel out, it does not have a large impact on volumetric change uncertainty, thus additional research on aleatory uncertainty is not a high priority. Increased understanding and better quantification of epistemic uncertainty, primarily from biases in depth soundings, could decrease the uncertainty in our estimates of change in sand storage in the bed of the Bay.