A Framework for Assessing the Bearing Capacity of Sandy Coastal Soils from Remotely Sensed Moisture Contents

Monthly field measurements were conducted from September 2020 to January 2021, in April 2021, in June/July 2021, and in October 2021 at the USACE Field Research Facility (FRF) in Duck, NC (Fig. 1) along the transect lines shown in Fig. 1(b). During this period, 35 individual surveys were conducted (Table 1) across all four seasons to determine if seasonal trends in bearing capacity, grain size characteristics, and moisture content existed. Sand samples and moisture probe measurements [Figs. 2(a and b), respectively] were conducted during all surveys, generally within a 1- to 1.5-h testing window. Satellite images were taken within the same testing window to avoid bias from changes in the tidal stage. Measurements along each transect stretched from the toe of the dune to the water line; due to fluctuations in water levels impacting the length of the emerged beach profile, the number of stations sampled per survey varied. Generally, stations were spaced 5 m cross-shore, and the three transects were spaced 20 m along-shore. Measurements on December 9, 2020 were spaced 7.5 m cross-shore, and measurements on November 13, 2020 were collected along a single along-shore transect [orange squares in Fig. 1(b)] due to high water conditions. During field measurements, heavy mineral materials were intermittently observed, especially in November–December 2020, and particularly around the second to fourth sampling stations. The presence of heavy minerals was easily identified from sand grain color and resulted in elevated unit weights and moisture content measurements from the moisture probe.

Historically, the beach near the USACE-FRF is composed of a predominately quartz sand that ranges from well-graded to uniform sand, with $d_{50}$ between 0.3 and 3.84 mm (Birkemeier et al. 1985; Gallagher et al. 2016; Stauble 1992). The site is typically a barred beach and morphological features such as foreshore beach cusps are known to form, typically following storms and during energetic fall, winter, and spring months (Holland 1998). In the present study, sediment sampling was limited by the restricted measurement window of 1–1.5 h. In total, 128 sand samples were carefully collected in 10-cm-long transparent tubes [Fig. 2(a)]; the number and location of samples collected per survey are presented in Table 1. Here, the lowest station number refers to samples near the dune, with the station number increasing while approaching the water line. Moisture contents, unit weights, and grain sizes were determined for each sample. Samples were divided into two 5-cm sections for drying, with the top 5 cm sieved following ASTM D6913-17 (ASTM 2017). At least two physical samples were collected during each survey, with the sampling locations presented in Table 1. Most samples were collected along Transect A, with other samples randomly collected along other transects when visible changes in grain size, mineralogy, or density occurred. At each station, 3–5 volumetric moisture content measurements were collected using a Dynamax SM150 moisture probe [Fig. 2(b)], with readings in mV. These were then converted to moisture content using custom calibration constants developed for the site-specific saltwater and quartz sand (Brilli et al. 2022). Calibration constants were not developed for the heavy mineral sand observed during Fall 2020.

The PFFP blueDrop [Fig. 2(c)] was deployed at 415 stations to measure the quasi-static bearing capacity (qsbc) during 18 surveys in September–December 2020, April 2021, and June/July 2021 (Table 1). Full details of the PFFP data processing can be found in the Supplemental Materials. The dynamic nature of the test indicates that a dimensionless strain rate correction factor needs to be applied; for sandy soils, the logarithmic function is commonly used for submerged sands (Stark et al. 2009, 2012; Stephan et al. 2015; Stoll et al. 2007) and will be applied here. For this study, $K$, a dimensionless coefficient used in strain rate correction, was varied between 1.0–1.5 and compared to computed values of bearing capacity of partially saturated soils. At each station, 3–5 deployments (i.e., individual, independent measurements) were conducted and spaced about 1 m apart, for a total of 2022 deployments as shown in Table 1, with 2004 analyzed. An example of deceleration profiles from the PFFP is shown in Fig. S1 for Station A2 on December 7, 2020. Individual qsbc profiles were calculated for each deployment, with the individual deployments at each station averaged together in 1-cm increments. Then, this average qsbc profile was used to represent each station, for a total of 404 stations. At nine stations, accompanying moisture probe measurements were not collected; these stations were excluded from further analysis. Stations where the penetration depth was $<7\mathrm{cm}$ were also excluded from further analysis (Albatal et al. 2020). At penetration depths $<7\mathrm{cm}$, uncertainty in qsbc increases due to assumptions related to the geometry of the conical tip and the shear failure planes; this resulted in the removal of 153 stations. Finally, stations where heavy minerals were present were also excluded, resulting in 23 stations being removed. These stations were selected based on visual observations in the field, inspection of physical samples, and when moisture probe measurements were significantly higher than the surrounding measurements. This yielded a total of 227 stations from 18 survey periods used for further analysis. In total, about 93.5% of these stations were found within the subaerial or intertidal zones, representing exposed sediments typically imaged by satellites.

For partially saturated soils, Vanapalli and Mohamed (2007) suggested that the soil suction ($u_a-u_w$) should be included in Eq. (2) in order to account for the effects of partial saturationwhere $S$ is the degree of saturation calculated from the porosity and measured moisture content; ${\varphi }^{\prime }$ = friction angle; and $\psi$ = bearing capacity fitting parameter (equal to 1 for sandy soils). For sands, $c^{\prime }=0$.

Grain size distribution curves of 127 of the samples collected over all surveys are shown in Fig. 5; see Table S1 for details. Samples were classified according to Wentworth (1922): fine ($d_{50}=0.125–0.25\mathrm{mm}$), medium ($d_{50}=0.25–0.5\mathrm{mm}$), coarse ($d_{50}=0.5–1.0\mathrm{mm}$), and very coarse ($d_{50}=1.0–2.0\mathrm{mm}$). In total, 36 (18.3%) samples classified as fine, 81 samples (63.5%) classified as medium, 9 (7.9%) classified as coarse, and one sample (0.8%) classified as very coarse sand. Following the USCS classification (ASTM 2017), the fine, medium, and coarse sands are poorly graded sands (SP) while the very coarse sample was a well-graded sand (SW). The mean ($\mu$), standard deviation ($\sigma$), and range of grain size characteristics for each classification are shown in Table 3. Since the very coarse sand class represented only one sample and was collected near the water line at a visible gravel patch at the time of measurements, it will not be discussed further. The cross-shore variation in $d_{60}$, $d_{50}$, and $d_{10}$ is shown in Figs. 6(a–c), respectively. The gray x’s represent individual measurements while the orange x’s represent the average values for each station. The variability in the grain size characteristics was found to be generally independent of the season and instead was more strongly related to the location of the beach where the measurements occurred. In general, the grain size characteristics remain approximately constant until Station A7. Samples from A7 onward were generally collected in the swash zone, where considerable variation in grain size has been previously observed (Stauble 1992). Overall, the measurements matched well with observations by Stauble (1992) that were collected monthly at Transect line 62 (modern-day FRF line 1006) are shown in Fig. 6(b) as stars.

The dry unit weight, porosity, and friction angles are shown in Table 3 for each sand class. In total, 105 samples were used to find the unit weight and porosity for each class; samples that exhibited heavy mineral materials and those without unit weight measurements were excluded. The porosity ranged between 0.404 and 0.448, which is within the expected range of 0.3–0.6 (Hillel 1998) and with those previously reported for beach sands (Atkins and McBride 1992). As the mean grain size increases, the porosity decreases and the dry unit weight increases as a result of smaller particles filling the pores between larger particles (Hillel 1998); these trends are reflected across the sand classes. The measured angle of repose falls within the expected range of 30°–45°, with 34° as a typical value for dry sand (Glover 1989) and increased with increasing grain size. Previous tilt table tests performed on samples from Duck indicated an angle of repose of 37° (Stark et al. 2017), which may be a result of the slightly coarser material observed and tested (minimum $d_{50}$ of 0.347 mm). Previous work and correlations (e.g., Decourt 1989; Duncan et al. 2014) indicated that the friction angle should increase with increasing grain size. The fine and medium sands exhibited the same values for the angle of repose, likely due to the similarities in the grain size distributions of the two classes (Fig. 5).

The results of the sensitivity analysis are shown in Fig. 7 as a function of the moisture content for the three sand classes and all sand data combined; the presented values represent the maximum percent differences. Overall, the friction angle (dashed line) had the most significant impact on the estimated bearing capacity, with percent differences of 32%–48.4% for the different sand classes and 81.0%–81.7% for all three sand classes together. This is a result of the friction angle appearing in multiple locations in Eq. (3); particularly, the bearing capacity factors [$N_c$ and $N_{\gamma }$; see Eqs. (S9) and (S10) for details] vary significantly as a result of the varying friction angle. For the total range of friction angles, $N_{\gamma }$ varies between 38.4 and 220.8 and $N_c$ varies from 42.2 to 118.4, which would result in a wide range of bearing capacities. For the friction angle, the percent difference did not exhibit a strong relationship with the moisture content, with the percent difference varying by $<1\%$ for all four cases. Only dry angles of response were used to represent the friction angle since this represents the minimum value. The angle of repose was not tested at a varying moisture content despite the known influence of moisture content (Stark et al. 2017).

The impact of the remaining four properties ($d_{60}$, $d_{10}$, porosity, and dry unit weight) on the estimated bearing capacity related noticeably to changes in the moisture content. Particularly, the $d_{60}$ had a pronounced relationship with the moisture content. The two parameters that are related to the $d_{60}$, $a$ [Eq. (5)] and $n$ [Eq. (6)], decrease with increasing $d_{60}$; $a$ and $n$ are then used to compute the suction stress [Eq. (4)]. The parameters $a$ and $n$ are not dependent upon the moisture content and only vary as a function of the $d_{60}$. The suction stress, however, is strongly related to the moisture content, thus impacting the bearing capacity. For small values of moisture content, the calculated suction stress (and, in turn, bearing capacity) increases for increasing values of $d_{60}$. At a moisture content of about 3.9%, however, higher suction stresses and bearing capacities occur for smaller values of $d_{60}$; this trend continues as the moisture content increases beyond 3.9%. Independent of $d_{60}$, suction stress decreases and the bearing capacity increases as the moisture content increases. Across all four sand cases, the percent difference associated with a varying $d_{60}$ (solid line with bars in Fig. 7) on the bearing capacity decreased while increasing the moisture contents to about 3.9%; at this moisture content, the smallest percent difference was observed (0.1%–0.4%). Then, the percent difference increased as the moisture content increased. When increasing moisture contents beyond about 3.9%, the calculated suction stress decreases with an increase in $d_{60}$. At a moisture content of 35%, the largest percent difference occurred for all sand classes combined (31.0%), followed by the medium (20.5%), coarse (12.8%), and fine (3.5%) and classes.

The impact of $d_{10}$ (solid line in Fig. 7), which is used to compute the parameter $n$ [Eq. (6)], also exhibited a pronounced relationship with moisture content, though to a lesser extent than the $d_{60}$. The parameter $n$, and in turn, the suction, decreases with an increasing $d_{10}$; this trend was consistent across the entire range of moisture contents. The largest percent difference (12.6%) was observed when considering data from all sand samples at a moisture content of 7.4%. For the fine and coarse sand classes, the maximum percent difference was 6.2% at moisture contents of 6.4% and 2.1%, respectively, while for the medium sand class the maximum percent difference was 2.5% at a moisture content of 7.1%. Overall, the $d_{10}$ had a lower impact on the bearing capacity than the $d_{60}$. The porosity was used as the saturated water content, ${\theta }_s$, in Eq. (4) (Fig. 4) to calculate the suction as well as the degree of saturation in Eq. (3). Across all four cases, variation in the porosity (dot-dash line in Fig. 7) led to increasing percent differences with increasing moisture content. For all four sand cases, the maximum percent difference when varying the porosity was less than 10%, indicating minimal impact. Finally, the bearing capacities estimated using a varying dry unit weight (dotted line in Fig. 7) exhibited percent differences that decreased with increasing moisture content. When varying the unit weight, the percent difference ranged from 13.3% to 19.0% for fine sand, 19.7% to 26.2% for medium sand, 19.1% to 20.7% for coarse sand, and 22.6% to 28.6% when considering data from all sand samples.

An example of the typical cross-shore trends in moisture content, peak qsbc at $K=1.25$, and penetration depth are shown in Fig. 8 for the survey on October 8, 2020 [Transects A, B, and C; see Fig. 1(b)]. Little variability in moisture content and bearing capacity was observed across seasons; instead, it was found that trends in bearing capacity and moisture content were strongly related to the cross-shore beach zone where measurements occurred and storm events. The moisture content [Fig. 8(a)] increased cross-shore toward the water line. The standard deviation between individual measurements at each station ranged from 0.05% to 16.6%, with an average of 1.7%, when excluding stations where heavy minerals were observed. Higher standard deviations were observed during storm events (1.8% to 16.6%, average of 6.9%) versus non-storm conditions (0.05% to 9.4%, average of 1.5%). High standard deviations were occasionally observed in the subaerial zone, especially from September to December 2020 when heavy mineral material was abundant and may have been included in the sand at the site at depth but was not observed on the surface or in the samples. The highest standard deviations in moisture content (average standard deviation of 4.2%) occurred at mean moisture contents of 10%–22.5%; these moisture contents generally occurred across the entire intertidal zone.

The qsbc [Fig. 8(c)] generally increased and penetration depths [Fig. 8(b)] decreased while approaching the water line, following trends observed previously (Brilli and Stark 2020; Reeve et al. 2018; Stark et al. 2022). Over all surveys, the impact velocity ranged between 2.31 and $5.93\mathrm{m}/\mathrm{s}$ and penetration depths ranged between 1.3 and 24.6 cm for individual deployments of the PFFP. Most surveys exhibited a range of penetration depths cross-shore, with deeper depths observed in the subaerial zone and shallower penetrations in the swash; examples of these trends from a single survey are shown in Fig. 8(b). Impact velocities varied somewhat between surveys, representing the larger user dependence of the device in air than in water. However, analysis resulted in similar magnitudes of qsbc between surveys despite differences in the impact velocities, and little variability in the impact velocity was observed (standard deviations of about $0.1–0.2\mathrm{m}/\mathrm{s}$ in most cases; 0.3%–27.8% of the mean, with an average of 3.6%) between deployments at a given station. Deployments in the swash exhibited the most variability in terms of PFFP measurements, likely a result of the rapidly changing conditions in this zone (e.g., grain size, unit weight, and swash uprush/backwash) occurring in this region (Dean and Dalrymple 2001; Stark et al. 2022).

Fig. 9 presents the moisture content measured using the moisture probe and the corresponding qsbc using $K=1.25$ for stations where sediment samples, moisture probe, and PFFP measurements were collected simultaneously (blue circles) and stations where only the moisture probe and PFFP were used simultaneously (orange squares). For stations where deployments of the PFFP and sediment sampling occurred simultaneously (29 of 61 stations of those included in analysis; moisture contents of 0.5%–15.8%), the values of qsbc computed using $K=1.0-1.5$ were compared to the estimated bearing capacity computed from the measured sand sample characteristics and the bearing capacity factors for cones from Durgunoglu and Mitchell (1973). The qsbc that resulted in the closest match was selected as the measured qsbc. The comparison between the measured and estimated bearing capacities are shown in Fig. 10(a), with a correlation coefficient, $\rho$, of 0.83.

Due to the excellent match between the measured qsbc and the bearing capacity predicted using measured sand characteristics, the bearing capacity was then calculated using assumed sand characteristics (227 of 404 total stations) and the measured moisture content at that station. For these stations, volumetric moisture contents ranged from 0.1% to 33.6%, the location of maximum qsbc ranged between 7 and 17 cm, and peak qsbc varied between 23.4 and 95.7 kPa. The medium sand characteristics measured in situ (Table 3) were used as the assumed sand characteristics to compute the bearing capacity for all stations. Medium sand represented 63.5% of all samples collected, and, when combined with fine sand (which exhibits similar characteristics to medium sand; Table 3) represented 81.7% of samples collected. From this, the measured qsbc and the bearing capacities estimated using the medium sand characteristics are shown in Fig. 10(b), with $\rho =0.71$.

To examine the model fit for estimating the bearing capacity using the medium sand class characteristics, the residuals of the estimated bearing capacities were computed, and are shown in Fig. 11. The residuals typically ranged between $-22.7\mathrm{kPa}$ and 33.1 kPa (negatives represent underestimates) with an average residual of 1.32 kPa, with 194 stations (or 85.5% of stations) exhibiting estimates within $\pm 10\mathrm{kPa}$. At one station, B7 on October 7, 2020, the bearing capacity was overestimated by 112.4 kPa (i.e., the bearing capacity model suggested stronger sands than measured by the PFFP) and exhibited the highest moisture content measured (33.6%). Previous studies found that PFFP measurements within the swash/water generally exhibit higher strengths and penetration depths $<7\mathrm{cm}$ (Brilli and Stark 2020; Reeve et al. 2018; Stark et al. 2022), typically a result of in situ densification from wave action (Heathershaw et al. 1981). However, station B7 exhibited a deeper penetration depth (15 cm) and lower measured strengths ($<50\mathrm{kPa}$) than typically observed within the swash. This reduction in strength may be a result of environmental factors such as the development of hydraulic gradients and swash dynamics, which have been observed to occur in sand beds under non-storm conditions (Horn et al. 1999; Stark et al. 2022). For sand with $d_{50}=0.24\mathrm{mm}$, Horn et al. (1999) found that large upward (negative) hydraulic gradients can occur in the uppermost layers (about 1–2 cm depth) of the bed under swash retreat, potentially leading to bed fluidization (Packwood and Peregrine 1980), reduced strengths, and increased PFFP penetration depths. This process was not considered in the model for estimating CI [Eq. (1)] nor in Fig. 4. To determine whether the development of hydraulic gradients could reduce the measured bearing capacity, a simple analysis was performed. For foundations subjected to seepage forces from groundwater flow, Veiskarami and Kumar (2012) suggested replacing the unit weight in Eq. (3) with ${\gamma }_{\mathit{sub}}\pm i{\gamma }_w$, where ${\gamma }_{\mathit{sub}}$ is the submerged unit weight and $i$ is the hydraulic gradient. Upward-directed flows decrease the bearing capacity (${\gamma }_{\mathit{sub}}-i{\gamma }_w$), while downward flows (${\gamma }_{\mathit{sub}}+i{\gamma }_w$) increase in the bearing capacity. Horn et al. (1999) found maximum upward gradients of 2.3 to 2.7 in the upper 1–2 cm of the sand bed, which are sufficient to fluidize the bed (Packwood and Peregrine 1980). These hydraulic gradients were tested in Eq. (3) by replacing $\gamma$ with ${\gamma }_{\mathit{sub}}\pm i{\gamma }_w$. This resulted in an estimated bearing capacity of 37.8 kPa, with a residual of $-5.5\mathrm{kPa}$ (slightly underestimated bearing capacity).

For the optical image, the estimated moisture contents are shown in Fig. 12(a) to the north of where field measurements occurred [Fig. 1(b)]. Lower moisture contents were observed near the toe of the dune and increased toward the water line, following previously observed trends (Paprocki et al. 2022). The estimated moisture contents for the SAR image are shown in Fig. 13(a) in the general vicinity of field measurements [Fig. 1(b)]. The SAR image confirmed the same trends in moisture content as the optic image. However, higher values of moisture content were observed near the toe of the dune due to the lower limit set by the SAR model (lower limit of $M_v=[-6.286/\ln (\theta /90){]}^{-1.538}$). The CI [Eq. (1)] was estimated using the estimated moisture contents from both the optical [Fig. 12(b)] and SAR [Fig. 13(b)] images. For the optical image, lower CI values were observed near the toe of the dune than for the SAR image, in line with the differences in the estimated moisture contents. In both images, higher CI values are observed near the water line where higher moisture contents are expected, consistent with the expected trends (Sullivan and Anderson 2000).

The estimated mean (solid black line) and standard deviation (dotted line) of the bearing capacity for the LARC-V is shown in Fig. 14 as a function of the moisture content, following the procedure shown in Fig. 4. Here, the parameters for medium sand from Table 3 were selected for calculating the bearing capacity since this class represented the majority (63.5%) of samples collected, and when combined with the fine samples (which exhibit similar characteristics to the medium samples, see Table 3), represented 81.7% of samples collected. It was assumed that the LARC-V was loaded to the gross weight of 8,618.3 kg and tires were inflated to a pressure of 68.9 kPa. Sinkage was calculated using the relationship shown in Fig. 15, which is based on field measurements of rut depth and moisture content during the October 2021 measurements. It should be noted that rut depth (and not sinkage) was measured, but it is assumed that they are equal since the difference between them is small (Affleck 2005). The applied bearing pressure was calculated as a function of the sinkage, and in turn, moisture content, ranging from 74.7–76.3 kPa, with lower values associated with lower moisture contents. The standard deviation of the moisture content was calculated using the equation derived from field measurements of the moisture probe [Eq. (S7)]. From Fig. 14, the estimated mean and standard deviation of the bearing capacity both increased as the moisture content increased.

Using the process shown in Fig. 4, the bearing capacity estimated from the moisture contents based on the optical image is shown in Fig. 12(c) and those based on moisture contents from the SAR image are shown in Fig. 13(c). For the optical image, where the gravimetric moisture contents were estimated, the void ratio was used for calculating the suction stress [Eq. (4)]. For the SAR image, porosity was used to calculate the suction stress since volumetric moisture contents were estimated. The most suitable paths for driving on the beach (based on the probability of inadequate bearing capacity) are presented in Fig. 12(d) for the optical image, and in Fig. 13(d) for the SAR image assuming a probability of failure of 10%. For the optical image [Fig. 12(d)], areas near the toe of the dune are unsuitable for trafficability, which corresponds to areas of lower bearing capacity [see Fig. 8(c)]. However, the SAR image [Fig. 13(d)] indicates that the entire area of the beach is suitable for driving, likely a result of the lower limit of the moisture content used when modeling the moisture content.

From sensitivity analysis (Fig. 7), the selected friction angle has the greatest impact on the bearing capacity since the friction angle explicitly appears in Eq. (3) and in the equations for the bearing capacity factors [see Eqs. (S9) and (S10)]. Therefore, care should be given toward selecting the sand class (since this dictates the friction angle) in order to not overestimate the strength. Overprediction of the friction angle would lead to potentially unsafe driving conditions and increase the risk of bearing capacity failure. Thus, if the literature suggests that both medium and coarse sand are present at the site, it would be conservative to select the medium sand characteristics. For example, Stauble (1992) [black stars in Fig. 6(b)] suggested that medium sand dominates the subaerial zone for non-storm conditions [a distance from the dune of 8 and 17 m in Fig. 6(b)], with slightly coarser material observed during storm events. This coarse material likely resulted from finer material being washed away during high water conditions during storm events and wind transport during fair conditions (Stauble 1992). Thus, the use of a lower friction angle in the subaerial would yield conservative estimates for bearing capacity. During non-storm conditions, the intertidal (black stars at about 30.6 m from the dune) exhibited significant variation in grain sizes, ranging from fine ($d_{50}=0.22\mathrm{mm}$) to very coarse ($d_{50}=3.81\mathrm{mm}$) sands. However, medium sand still dominates [about 45.5% of all samples observed in Stauble (1992)], indicating that it is still reasonable to use the medium sand characteristics to predict the strength in the intertidal zone. From Glover (1989), the angle of repose for loose sands should range between 30° and 45°, with 34° a typical value for dry sands. Therefore, in the absence of grain size information, friction angles of 30°–34° can be used to provide a conservative estimate of the bearing capacity.

When using the bearing capacities derived from assumed sand characteristics using data from in situ measurements [characteristics for the medium sand class from Table 3; Fig. 10(b)], there is general agreement between the measured qsbc and estimated bearing capacities, suggesting that the presented approach (flow chart in Fig. 4) can be applied for average (i.e., non-storm) beach conditions. While 194 stations had residuals smaller than 10 kPa (Fig. 11), the remaining 33 exhibited residuals outside this range. Of these, 11 exhibited negative residuals, indicating that the bearing capacity was underestimated, representing conservative results. The remaining 22 exhibited overestimates of the bearing capacity to some degree. In total, 14 of these stations were located within the subaerial zone during normal (non-storm) conditions. The subaerial zone is known to exhibit low relative densities and unit weights, with unit weights likely lower than the average values used here (Table 3). NAVFAC (1986) suggests that loose, uniform, fine-to-medium sands, which would typically characterize the subaerial zone, can have unit weights as low as ${13.2\mathrm{kN}/\mathrm{m}}^3$. In this study, 38% of medium sand samples collected at Station A1 exhibited moist unit weights less than this when using the average dry unit weight to estimate the moist unit weight. This indicates that it may be conservative to use slightly lower mean values of dry unit weight in this region in order to not overestimate the strength, especially given the influence of the unit weight on the bearing capacity (Fig. 7). Additionally, four measurements with residuals $>10\mathrm{kPa}$ occurred in December 2020 when a significant amount of heavy mineral material was observed near the toe of the dune. Two of those measurements exhibited moisture contents (8.0%–11.3%) above the expected range (Paprocki et al. 2022), pointing at the presence of heavy mineral material within the soil that was not visible, resulting in elevated moisture contents and overestimated bearing capacities. Finally, organic materials were observed at Station B4 on September 9, 2020, likely elevating the measured moisture content.

High residuals were also observed at stations where moisture contents vary significantly. One such station was located at the transition between the subaerial and intertidal zones (where the moisture content can vary rapidly over short distances), while the other two were measured near the edge of the water in the lower intertidal zone (where the moisture content at the edge of the water can vary rapidly in response to wave runup and tides). Additionally, four of those measurements were collected during storm conditions on November 12, 2020 [about 7.5–8.5 cm of rainfall over 48 h (NOAA 2021b)]. For these six stations (four during storm, two near the edge of the water), wetting occurred during measurements, and the framework for estimating the bearing capacity may not be applicable to these measurements; Eqs. (5) and (6) were developed for drying (and not wetting) conditions to estimate the suction stress. During in situ measurements, PFFP deployments and moisture probe measurements were taken within minutes of one another, but even during this short time period large variations in moisture contents were observed. This highlights the dependence of the bearing capacity on moisture content as demonstrated by the sensitivity analysis (Fig. 7), suggesting that accurate moisture contents and consideration of variability are needed to safely estimate the bearing capacity.

Another factor likely influencing the measured bearing capacities within the intertidal and swash zones is the variation of grain sizes. Particularly, the $d_{60}$ had a significant impact on the estimated bearing capacity in the sensitivity analysis through its relationship to the moisture content (Fig. 7). At higher moisture contents ($>20\%$), the bearing capacity was almost always overestimated (Fig. 11), corresponding to when a larger percent difference was observed in the medium sand class [Fig. 7(b)]. Throughout the survey sequence, surficial gravel (coarser sediments) as well as a wider range of grain sizes (Fig. 6) were found near the water line (within the intertidal and swash zones). This trend has been observed in previous studies of sandy beaches (Abuodha 2003; Dingler and Reiss 2002; Gallagher et al. 2016, 2011; MacMahan et al. 2005; Prodger et al. 2016, 2017; Reniers et al. 2013; Stark et al. 2022; Stauble 1992; Thornton et al. 2007). Even for beaches that are usually composed of fine and medium sands, larger grains have been observed during and following storm events (Dingler and Reiss 2002; Gallagher et al. 2011; Prodger et al. 2016, 2017; Stauble 1992), indicating that the characteristics of coarse sands from Table 3 may be applicable at the water line or during storm events. Future work will focus on spatially mapping the distribution of grain sizes within the swash zone to improve bearing capacity estimates in this region.

The remaining station that exhibited a residual $>10\mathrm{kPa}$ was B7 on October 7, 2020, which was located in the swash at the time of measurements. While larger grain sizes than those found for the medium sand class likely occurred in this region (Fig. 6), when a $d_{60}$ of 3.81 mm [the largest $d_{50}$ observed by Stauble (1992) in the intertidal zone] was tested, the bearing capacity was still overestimated. This indicates that, despite the importance of the grain size on the estimated bearing capacity, other external forces, such as upward-directed hydraulic gradients and swash processes, are likely key processes influencing the qsbc near the waterline. The impact of these pressure gradients can be captured with a simple calculation using values presented in the literature to represent a theoretical “worst case” scenario for bearing capacity within the swash. Of note, the large negative gradients measured by Horn et al. (1999) used in this study (2.3–2.7) occurred during flood (rising) tide conditions, with positive (downward) gradients dominating during ebb (falling) tides, indicating that the tidal stage also has an impact on the development of pressure gradients that may reduce the bearing capacity. Measurements in October 2020 began at low tide with Transect A, with the final transects, B and C where the low strengths were observed, occurring during the start of flood tide (NOAA 2021a). Thus, it may be suggested to use upward-directed hydraulic gradients when estimating the bearing capacity as a conservative boundary during flood tide. In the future, measurements of the bearing capacity and pressure gradients should be collected simultaneously to better understand their relationship.

Uncertainty in estimating the bearing capacity may arise due to three sources: inherent variability of the soil properties due to natural processes altering the soil, measurement error, and transformation uncertainty (Phoon and Kulhawy 1999). To account for the inherent variability of soil properties, a probability of inadequate bearing capacity [Eqs. (10)–(12)] was computed, assuming a normal distribution. Baecher and Christian (2005) found this to be a reasonable assumption, likely resulting in an overestimation of the probability of failure. In Eqs. (10)–(12), it was assumed that the input parameters followed normal distributions. This is a common assumption for geotechnical properties, as most geologic processes follow either a normal or lognormal distribution (Kadar and Nagy 2017; Lacasse and Nadim 1998). Lacasse and Nadim (1998) reported that the friction angle is normally distributed, with coefficients of variation (COV, $\sigma /\mu$) of 2%–5%. Lacasse and Nadim (1998) also reported that, for all soils considered together, the submerged unit weights and porosities were also normally distributed, with COVs of 0%–10% and 7%–30%, respectively. Duncan (2000) found COVs of 3%–7% for the unit weight and 2%–13% for the effective stress friction angle when analyzing the data reported from a number of sources. Finally, From Table 3, the COVs range from 8.9% to 12.5%, 6.1% to 8.8%, and 2.1% to 2.8% for the porosity, unit weight, and friction angle, respectively, and are within the expected ranges from Lacasse and Nadim (1998) and Duncan (2000). Phoon et al. (1995) suggested COVs of 5%–11% for the effective stress friction angles of sands, and are higher than those presented in Table 3; however, the test method was not specified for these data, and thus, a direct comparison cannot be made. Limited data is available regarding the distribution of grain sizes, though Kadar and Nagy (2017) assumed normal distributions for $d_{60}$ and calculated a COV of 26% for the $d_{60}$ of a sand/silty sand; this value similar to the values observed for medium (29.4%) and coarse (39.1%) sands. However, further testing is required to investigate the distributions of the grain sizes at sandy beaches. Lacasse and Nadim (1998) found that varying the distribution function of soil parameters had a minimal impact on the computing the probability of failure, indicating that the selection of normal distributions for sediment characteristics is an acceptable assumption.

For both images, moisture contents were averaged over a $3\times 3\mathrm{pixel}$ area, or an area of $1.5\times 1.5\mathrm{m}$ for the SAR image and an area of $3.72\times 3.72\mathrm{m}$ for the optical image. During field measurements, each station had five measurements collected that were spaced about 1 m apart. Most stations (for non-storm measurements and when stations with heavy minerals were excluded, 523 of 609 stations) exhibited limit variability (standard deviations $\le 1.5\%$) between the five measurements of moisture content, indicating that the resolutions of the images used are likely suitable to estimate the moisture content and in turn, sediment strengths. Additional scrutiny should be given to estimates at the edges of beach zones and in the presence of geomorphological features, where more variation was observed between the individual measurements at each station.

The estimated CI values [Figs. 12(b) and 13(b) for the optical and SAR images, respectively] indicate that the entire area of the beach would be suitable for trafficability, as a minimum CI of 60 is required for the LARC-V to safely traverse sandy beaches (Jenkins 1985). However, this approach does not account for variable soil properties (i.e., changes in grain sizes, unit weights, friction angles, and porosities) beyond the moisture content. The estimated probability of failure [Figs. 12(d) and 13(d)], which does incorporate this variability, can be used to select the best path for driving across the beach and contribute to risk reduction. For both the high-resolution optical and SAR data, both the estimated moisture contents and bearing capacities followed the expected trends from field measurements (Fig. 8), with lower strengths observed near the toe of the dune and higher strengths near the water line. The estimated moisture contents from the optical image [Fig. 12(a)] spanned the entire range of possible moisture contents, with lows of about 0% and peaks of 25%. Near the toe of the dune, there were high estimates of moisture content (about 10–15%), which may be a result of both increased surface roughness (Dean and Dalrymple 2001) and organic material, both which are known to decrease the reflectance and impact the estimated moisture contents (Baumgardner et al. 1986; Johannsen 1969). This results in potentially erroneous estimates of bearing capacity, and in turn, maps of safe driving. Thus, it may be reasonable to use a lower mean moisture content (i.e., about 2–3 standard deviations below the mean) for this region in order to more accurately assess the best path for driving in this region.

For the SAR image, both the estimated moisture contents [Fig. 13(a)] and bearing capacities [Fig. 13(c)] near the toe of the dune were higher than those estimated from an optical image. For the model used here (Paprocki 2022b), there is a lower limit for the estimated moisture contents that is a function of the incidence angle of the image (${M_v=[-6.286/\ln (\theta /90)]}^{-1.538}$ ; (Oh 2004)). This lower limit is 6.8% (based on $\theta =30.57°$ ) for the image used in this study, which would result in a bearing capacity of 119.2 kPa and a probability of failure of 2.1% (Fig. 14). The minimum moisture content measurement from October 8, 2020 [Fig. 8(a)] was 0.8%, which results in a bearing capacity of 76.7 kPa and a probability of failure of 46.8% (Fig. 14). From the October 8, 2020 measurements, about 48% of all measured moisture contents fall below the lower limit of 6.8%. These moisture contents occurred in the subaerial zone, which is generally characterized by lower strengths and unit weights. Therefore, in order to not overestimate the strengths in this zone, an alternative approach for predicting the most suitable path for driving on the beach is shown in Fig. 16 for the SAR image. Here, regions where the estimated moisture content is $<6.8\%$ are classified as “unsafe for driving” and regions where the estimated moisture contents were $>6.8\%$ are classified as “safe to drive.” While this map is not an explicit representation of strength, it is guided by a need to not overpredict the suitability of soils for trafficability and instead address the most suitable path for driving based on estimated moisture contents.

The two images selected for estimating the moisture content and strengths are representative cases of the range of images tested. For the optical image from March 14, 2018, the estimated moisture contents fell well within the expected range. Similar trends were observed for all other images tested, including different grain sizes, resolutions, and sites (Paprocki et al. 2022), indicating that the approach proposed in Fig. 4 can be used to estimate the bearing capacity for other high-resolution optical images with NIR bands. During measurements in 2020–2021, 28 X-band SAR images were collected, with a further five collected in 2018–2019, including three collected at approximately the same time as the optical images. For the model presented (developed for $\theta =30°-46°$, see the supplemental materials for details), similar trends of moisture contents as seen in Fig. 13(a) have been observed across a range of seasons, indicating that the presented approach can likely be applied for a range of seasons.

As a next step, the use of other satellite sensors to estimate the moisture content should be tested. In this study, only high-resolution X-band and multispectral images were used. C-band (wavelengths of 3.5–7.5 cm) data from satellites such as Radarsat-2 (ground resolutions of about $2\mathrm{m}/\mathrm{pixel}$) are an alternative candidate for estimating the moisture content, especially due to the similar ground resolutions to optical images tested in Paprocki et al. (2022). Lower resolution satellites such as the C-band Sentinel-1 SAR and the optical Sentinel-2 sensors (resolutions of about $10\mathrm{m}/\mathrm{pixel}$) likely do not have the appropriate resolutions to detect the subtle changes in moisture content that occur at sandy beaches of similar widths (on the order of tens of meters depending on tidal stage), especially given the dependence of the bearing capacity on accurate measurements of the moisture content (Figs. 7 and 11). However, they may be able to provide information regarding changes in ground conditions at beaches, such as grain sizes, beach slopes, and water lines. Sentinel-1 and -2 are publicly available and have revisit times of 5 days, meaning that data related to ground changes can be collected regularly to inform the detailed moisture content analysis from the high-resolution X-band SAR and WorldView optical data. Particularly, Sentinel-2 data has been used to map shoreline positions (Vos et al. 2019a, b) and beach slopes (Vos et al. 2020), which are important factors toward identifying the location of bare soils and estimating the moisture contents and, in turn, strength.