Earthquake Hazard Uncertainties Improved Using Precariously Balanced Rocks

Probabilistic seismic hazard analysis (PSHA) is the state‐of‐the‐art method to estimate ground motions exceeded by large, infrequent, and potentially damaging earthquakes; however, a fundamental problem is the lack of an accepted method for both quantitatively validating and refining the hazard estimates using empirical geological data. In this study, to reduce uncertainties in such hazard estimates, we present a new method that uses empirical data from precariously balanced rocks (PBRs) in coastal Central California. We calculate the probability of toppling of each PBR at defined ground‐motion levels and determine the age at which the PBRs obtained their current fragile geometries using a novel implementation of cosmogenic 10Be exposure dating. By eliminating the PSHA estimates inconsistent with at least a 5% probability of PBR survival, the mean ground‐motion estimate corresponding to the hazard level of 10−4 yr−1 (10,000 yr mean return period) is significantly reduced by 27%, and the range of estimated 5th–95th fractile ground motions is reduced by 49%. Such significant reductions in uncertainties make it possible to more reliably assess the safety and security of critical infrastructure in earthquake‐prone regions worldwide.


Introduction
Earthquakes can be profoundly destructive to societal infrastructure, human life, and the environment. It is estimated that over the past 40 yr the global population exposed to a moderate to severe intensity earthquake has increased 93% to 2.7 billion people (Pesaresi et al., 2017). The standard method used to assess the hazard for such potentially damaging earthquakes is probabilistic seismic hazard analysis (PSHA), which is a framework that allows the estimation of the rate or probability of exceeding a level of ground-motion intensity at a site (Gerstenberger et al., 2020, and references therein). What PSHA currently lacks is a quantitative method to use empirical data to both validate and refine the estimates of ground-motion intensities for long return periods used for seismic design that exceed the timescales of historical observations of ground shaking or damage to buildings. Geologic constraints are the only empirical data available to validate PSHA estimates over prehistoric timescales because the length of the earthquake cycle precludes the use of historical observations to assess longer-term patterns of seismicity. The development of robust and quantitative validation and evaluation methods to reduce uncertainties in earthquake ground-motion estimates are particularly required at long return periods of 10 3 to 10 6 yr because PSHA estimates for such return periods are highly uncertain and yet essential for the siting, design, and continued maintenance and monitoring of critical facilities, such as large dams, nuclear power plants, and nuclear waste repositories.
The validity of PSHA model estimates relies on the combination of accurate input parameters, including the characterization of proximal earthquake sources, that is, seismic source models, as well as predictive models of ground shaking, that is, ground-motion models (GMMs) such as NGA-West2 (Anderson & Brune, 1999;Bozorgnia et al., 2014;SSHAC, 2012). The various alternative input parameters and alternative models are incorporated within a logic-tree framework in which branches are weighted based on available data, physical constraints, and expert interpretation (SSHAC, 2012). The hazard is then calculated for each combination of the branches from different logic-tree nodes. To incorporate all the uncertainties in the PSHA model, some modern PSHA studies have over 1 million end branches of the logic tree (Stepp et al., 2001), with each end-branch combination producing its own hazard curve. The hazard curves from all the logic-tree branches are used to compute the expected (mean) hazard curve, as well as hazard fractiles, for example, 5th, 16th, 50th, 84th, and 95th, to capture the scientific (epistemic) uncertainties in the hazard estimates (SSHAC, 2012). At low rates of exceedance, the impact of differences between the alternative models becomes more significant because the computed hazard at these rare levels is sensitive to the modeling of the tails of the input parameter distributions.
These large uncertainties at low rates of exceedance have particular importance in terms of both ensuring adequate safety and potential design costs for critical infrastructure and facilities globally. For example, in the PSHA model for the previously proposed U.S. high-level nuclear waste repository at Yucca Mountain, the 5th-95th fractile peak ground accelerations (PGAs; in units of gravity, g) corresponding to the 10 −4 yr −1 hazard level range from 0.2 to 0.9 g (Stepp et al., 2001). The recurrence intervals of rare, large earthquakes generally exceed the observation length of instrumental records, so such records do not provide complete information on the surrounding seismic sources (Anderson et al., 2011). Analytical modeling of the physical limits of earthquake ground motions has previously been investigated as an approach to place constraints on large, infrequent earthquakes (Andrews et al., 2007); however, the ground-motion values constrained by these physical limits are too large to be useful for most engineered structures. Until now, it has not been possible to provide a quantitative validation of each of the individual output hazard curves to directly revise the PSHA model results and reduce the uncertainties of the PSHA estimates for the hazard levels and ground-motion values of greatest interest.
To validate such PSHA model estimates at low rates of exceedance, ground-motion constraints over long timescales are required in proximity to major earthquake sources. Ancient fragile geologic features have been previously identified as potentially useful for validating unexceeded ground motions estimated from PSHA models both in the United States (Anderson et al., 2011, and references therein;Hanks et al., 2013) and New Zealand (Stirling & Anooshehpoor, 2006). The development of a probabilistic method to analyze the unexceeded ground motions represented by such features was a significant advancement in the approach Hanks et al., 2013;. These methods, however, did not address how to apply the constraints from such features to validate each of the individual output hazard curves, identify the inconsistent parts of the PSHA model, and, in turn, directly revise the PSHA model and reduce the uncertainties in the hazard and resulting ground-motion estimates. There are several ongoing efforts to develop methods to test and update existing PSHA models, including using data from fragile geologic features (Stirling et al., 2019), but a method has never been publicly available.
To provide the data to use to update the PSHA results, we use a class of fragile geologic features known as precariously balanced rocks (PBRs) (Brune, 1996;Brune & Whitney, 2000), which are boulders balanced on, and mechanically separated from, a sub-horizontal pedestal and possess a susceptibility to toppling when exposed to earthquake ground shaking. For constraining hazard using PBRs, we need to estimate the likelihood of a PBR toppling due to a strength of earthquake shaking and the length of time that the PBR has been in its current fragile configuration. We use existing models to estimate the probability of failure, that is, toppling or overturning, given an intensity of ground shaking of each PBR, which we define as the fragility of the PBR (Anooshehpoor et al., 2004;Purvance, Anooshehpoor, et al., 2008). The time elapsed since each PBR achieved its current shape, and thus fragility, defines the fragility age of the PBR. In previous studies, a key limitation of using PBRs to constrain hazard has been the large uncertainty in the fragility age. We use a new bespoke implementation of cosmogenic isotope ( 10 Be) surface exposure dating (cf. Balco et al., 2011) to develop improved PBR fragility ages. Following the methodology developed for Yucca Mountain Hanks et al., 2013), we use the PBR fragility and fragility age in combination to provide upper bounds on past ground motions that occurred at a site in coastal Central California by either onshore or offshore seismic sources. Our study site of Double Rock (Figure 1a) is located close to both the population center of San Luis Obispo and the critical infrastructure of the Diablo Canyon nuclear Power Plant (DCPP). With our improved method for estimating the fragility age, we develop the most rigorous method to date of using PBR constraints to reduce uncertainties in the hazard estimates and thus improve the PSHA model at long return periods.
This new method of validating all the hazard outputs of a PSHA model and reducing uncertainties in PSHA estimates is an important step forward in broadening the use of PBRs to improve the reliability of PSHA models. The validation of PSHA estimates using the observed geologic PBR data enables the elimination of combinations of seismic source model and GMM alternatives within the PSHA model that are inconsistent with the age and fragility of the preserved PBRs. The reduction in uncertainties provided by our analysis of PBRs delivers significantly more precise seismic hazard assessments, particularly over timescales greater than 10,000 yr for which uncertainties in hazard estimates are greatest. The method we present in this study can be applied to all sites where PBRs are preserved, including, and most importantly, where large uncertainties in hazard assessments and critical infrastructure or human populations coexist .

Tectonic, Geologic, and Geomorphic Context
The seismic hazard at our study site is dictated by the tectonic setting of the region. The seismotectonics of coastal Central California result from the relative motion between the Pacific and North American plates, which is estimated at 50.2 ± 1.1 mm yr −1 toward N35.8°W ± 0.2 (1σ) by the MORVEL plate motion model (DeMets et al., 2010). In the study area, plate motion is oblique to the plate boundary and partitioned between dextral strike-slip faulting and reverse and oblique-reverse-slip faulting west of the San Andreas fault (Clark et al., 1994). At the Double Rock site in coastal Central California, the seismic source contributing to the majority of the earthquake hazard is the Hosgri fault, which is an offshore dextral strike-slip fault 6 km southwest of the site (Figure 1a). The Hosgri fault is part of the plate boundary fault system and has a slip rate of 0.6-3.1 mm yr −1 with a mean characteristic magnitude of 7.1 (Pacific Gas and Electric Company [PG&E], 2015b). Other proximal, less hazardous faults included in the Double Rock site PSHA model are the San Luis Bay fault, with a slip rate of 0.09-0.34 mm yr −1 , the Los Osos fault, with a slip rate of 0.13-0.56 mm yr −1 , and the Shoreline fault, with a slip rate of 0.03-0.16 mm yr −1 (PG&E, 2015b). All these seismic sources are near to our study site ( Figure 1a) at distances that could lead to significant ground shaking during large magnitude earthquakes. There are epistemic uncertainties, however, in the characterization of these seismic sources and their ground motions. For example, in the logic-tree input to the PSHA model for the DCPP, the Hosgri fault is given branches for alternative slip rates of 0.6, 1, 1.7, 2.3, and 3.1 mm yr −1 . There are also branches for different models of the amplitude of shaking for a given earthquake scenario: For example, the median peak ground acceleration from a M7 earthquake at a distance of 5 km used for DCPP ranges from 0.27 to 0.50 g (GeoPentech, 2015).
In order to test the PSHA model (Abrahamson et al., 2019;PG&E, 2015a) estimates for the location of the Double Rock site, we utilize PBRs preserved on coastal landforms. Tectonically uplifted marine terraces are commonly preserved along the coast of California and have been used to estimate a Pleistocene (post-120 ka) uplift rate of~0.2 m kyr −1 in the region of the Double Rock site . Double Rock East (DRE) and Double Rock West (DRW) are abandoned sea stacks that extend above two tectonically uplifted coastal marine terrace platforms that are separated by a paleo-shoreline ( Figure 1b). The platforms on either side of this paleo-shoreline are dated to~80 and~125 ka, which correlate to previous interglacial sea-level highstands during marine oxygen-isotope sub-stages 5a and 5e, respectively . DRE and DRW are composed of Franciscan Complex chert bedrock (Dibblee & Minch, 2006); the high resistance to weathering and erosion of chert relative to other rock types at the site results in the formation of sea stacks that extend above their respective highstand marine platforms. The chert bedrock is also characterized by near-vertical bedding, joints, and extensive deformation fabrics. These fabrics produce meter-scale blocks, some of which are detached from the sea stack outcrops to form PBRs. The processes of PBR formation are observed acting on a modern offshore chert sea stack (DRM, Figure 1b and Supporting Information Figure S1), which is within the wave-action zone of the current interglacial sea-level highstand. Meter-scale chert blocks eroded from the modern sea stack and deposited onto the modern platform indicate that the sea  Dibblee & Minch, 2006) overlain on LiDAR hillshade image. The three chert outcrops on which the fragile geologic features form are labeled, and the localities investigated in this study are numbered. Field photo inset (looking toward 140°) shows the paleo-sea stacks DRE and DRW extending above the abandoned major alluvial fan complex Qf1 and modern sea stack DRM.

10.1029/2020AV000182
stacks in the wave-action zone rapidly erode by block removal processes. Thus, DRM provides an analog for when the DRE and DRW paleo-sea stacks were being actively eroded during the 5e and 5a highstands.
Relative sea-level fall following the 5e and 5a highstands, in combination with tectonic uplift, resulted in platform abandonment and preservation of the DRE and DRW paleo-sea stacks . Stability of the landscape following platform abandonment, and therefore sea stack abandonment, has allowed for the preservation of any PBRs that formed by post-abandonment weathering and erosion processes. In response to platform abandonment and associated sea-level fall at the end of each highstand, thick (20-30 m) alluvial fan sediments sourced from now inactive large catchments to the northeast were deposited on the abandoned marine terraces . Our observations at the site confirm subsequent inactivity and thus stability of the major alluvial fan complex (Figure 1b) as demonstrated by channel incision both to the north and south of the Double Rock chert paleo-sea stacks, which suggests long-term preservation of any PBRs that have formed since major fan abandonment. Following the abandonment of the major fan complex, a set of smaller minor alluvial fans (Figure 1b) sourced from small catchments in the Irish Hills to the northeast were deposited on top of the older major fans; however, these minor fans did not extend to DRE or DRW. Colluvium deposits have also subsequently accumulated in small volumes on the abandoned major alluvial fan surface in the areas surrounding the paleo-sea stacks; this colluvium volume is predominantly meter-scale blocks that indicate erosion of the chert outcrops by block removal to form PBRs after the time of major fan abandonment.

Fragility
Field studies were conducted at the Double Rock site to comprehensively survey the chert paleo-sea stack outcrops, DRE and DRW, for candidate PBRs (Figure 1b). The seven best PBRs at the site were selected for detailed analysis based on the suitability of geomorphic development by block removal processes, block geometries that allowed for accurate fragility characterization, and their accessibility on the paleo-sea stacks (section 5). The geometric parameters that are required for our fragility characterization method were derived from our field observations and analyses of the geometry of these seven PBRs (Figures 2a and S2a  and Table S1).
Earlier methods of seismic hazard evaluation using PBRs applied a deterministic approach, in which the PBR fragility was defined as a single threshold ground motion required for toppling. Under this previous definition, the calculated quasi-static (Weichert, 1994) and dynamic toppling accelerations (Anooshehpoor et al., 2004) are estimates of the ground motions that have not been exceeded since the PBRs became fragile. PBRs with toppling accelerations of 0.1-0.3 g were classified as precarious, and those with toppling accelerations 0.3-0.5 g were classified as semi-precarious (Brune, 1996). The single deterministic PGA fragility assigned to a PBR was directly compared to hazard maps and curves (Anderson & Brune, 1999;Anderson et al., 2011). Using these previous methods, the single deterministic dynamic thresholds (acceleration that causes toppling) of the studied PBRs at the Double Rock site range from a minimum of 0.23 g for the most fragile PBR, DRW1 (Figure 2a), to a maximum of 0.64 g for the least fragile PBR, DRW5. We present these deterministic thresholds to allow direct comparison of the degree of fragility of the seven Double Rock PBRs to previous PBR studies. These thresholds, however, are not best suited for validating PSHA hazard curves because not only do they not consider the probability of the survival of the PBR at a given level of shaking, but also they do not evaluate the epistemic uncertainty in the estimation of the probability of survival.
To address this uncertainty in the survival of the PBR, we follow the Hanks et al. (2013) and Baker et al. (2013) methodology and define the fragility of a PBR as the probability of toppling given a strength of ground shaking (Purvance, Anooshehpoor, et al., 2008); since the PBR remains untoppled, there is a high probability that no previous earthquake shaking has exceeded this strength of shaking in the time that the PBR has been in its present fragile state. The probability of toppling a PBR during an earthquake event is a function of the PBR's calculated geometric factors (Table S1), the amplitude of the ground-motion excitation force, and the associated waveform of the force. We employ an existing fragility model based on the normalized geometric shape of a PBR that provides the toppling probability as a function of PGA and ratio of peak ground velocity (PGV) to PGA (PGV/PGA) (Purvance, Anooshehpoor, et al., 2008;section 5.1.3). This fragility model assumes that a PBR has a 2-D geometry, oscillates about two rocking points, remains infinitely rigid, impacts inelastically, does not slide, and conserves angular momentum on impact (Housner, 1963). This model has been validated by shake table experiments (Purvance, Anooshehpoor, et al., 2008), which show that a 2-D geometry generally underestimates the rocking response compared to that of its 3-D geometry (Anooshehpoor et al., 2004;Purvance, Anooshehpoor, et al., 2008); therefore, the 2-D geometry is a conservative estimate of the PBR's fragility, and the 3-D geometry may be more fragile. This fragility model also describes the rocking response of PBRs with geometric factors that are both symmetric and asymmetric (Purvance, Anooshehpoor, et al., 2008). We classify the PBRs at Double Rock as symmetric when the geometric factors that control the probability of toppling are within 5% of the values calculated for both fragile directions (section 5.1.2; Table S1).
The probability of toppling of each of the seven PBRs in this study was calculated across a range of values of both PGA and PGV/PGA. The two ground-motion parameters are used to represent both the amplitude of the short-period shaking (PGA) and the amplitude of long-period shaking relative to the short-period shaking (PGV/PGA). Purvance, Anooshehpoor, et al. (2008) showed that the predominant pulse duration of PGV/PGA is a strong indicator of overturning because of the increased period of time that a PBR is forced to tip in one direction. The probability of toppling is shown for the most fragile PBR at the Double Rock site, DRW1 (Figure 2b), which shows that in this fragility model, increased values of PGA and PGV/PGA both lead to an increased probability of toppling. Fragility results for the other six PBRs are provided in Figure S2b.

Fragility Ages
To estimate the fragility ages of the PBRs using our bespoke implementation of cosmogenic isotope ( 10 Be) surface exposure dating (Balco et al., 2011), a geomorphic model for PBR development at the study site is required. The PBRs, and their fragility, at the Double Rock site are primarily the result of erosion of the rock outcrop by meter-scale block removal, in which blocks of rock become detached through weathering and mass wasting processes. From our field observations, we interpret that the PBRs have developed as blocks formed between the joints, the steeply dipping bedding planes, and the deformation fabrics within the outer layers of the chert sea stack rock outcrops. At the outcrops we studied, the fabric in the chert bedrock creates blocks with an average height of 2.56 m, width of 2.13 m, and depth of 0.68 m ( Figure S3). Some of these blocks remain in place as PBRs after more unstable portions of the rock outcrop have fallen away around them. For a theoretical sample point located behind and shielded by a large block of average dimensions, we estimate that 92.7% of the incoming cosmic ray flux will be shielded by the block and surrounding rock outcrop ( Figure S3). Therefore, we assume the removal by erosion of such blocks from the rock outcrop exposes new surfaces that contain insignificant cosmogenic nuclides at our sample points (sections 3 and 5.2.1). Consequently, sampling of such rock surfaces on all sides of the PBRs allows us to date the timing of block removal events responsible for establishing the current PBR fragility and, in turn, determine bounding ages for the development of its current fragility (Figures 2c and S4). Our geomorphic model also assumes negligible cosmogenic 10 Be inheritance (nuclides produced prior to sea stack abandonment) and negligible erosion of accumulated 10 Be after surface exposure by block removal. Since both assumptions could affect our interpretation of exposure ages from 10 Be concentrations in sampled rock surfaces, we conduct analyses to test these assumptions by directly quantifying this inheritance and erosion.
First, to verify our assumption that the cosmogenic nuclide concentration at the time of platform abandonment (i.e., inheritance) is negligible, we analyzed three chert rock surface samples collected from the modern wave-cut notch in the active intertidal zone ( Figure S1). The average 10 Be concentration in the three modern samples is 5,280 atoms g −1 , with a standard deviation of 919 atoms g −1 (Table S3). This concentration amounts, on average, to 7% of the 10 Be concentrations in the paleo-sea stack block removal exposure age samples. With these data, we demonstrate no significant inherited 10 Be concentration in the chert prior to sea stack abandonment; nevertheless, all reported ages have been corrected for this inheritance value and uncertainties propagated in quadrature. Due to this inheritance correction, only the 10 Be inventory accumulated after sea stack abandonment is used in our exposure age calculations.
Second, to quantify the effect of our assumption of negligible erosion after surface exposure by block removal on our interpretation of exposure ages, we calculate the average erosion rate of the sea stack rock outcrops, and likewise an estimate of the rate of PBR erosion, from the volume of chert colluvium accumulated on the major alluvial fan surface since the time of fan abandonment (see sections 3 and 5.3.2). The volume of colluvium on the major alluvial fan surface where the PBRs are preserved on the southern side of the DRE and DRW outcrops is~67 and~150 m 3 , respectively ( Figure 1b). The minimum time since the abandonment of the alluvial fan surface is constrained by a fallen block (DRB-2) dated at 49 ± 4 ky (1σ uncertainties throughout) ( Figure S5), which is consistent with observations of the soil profile development (section 5.3.1 and Text S1). Using these volumes and age together gives an average maximum erosion rate of the paleo-sea stacks DRE and DRW of 2.4 and 2.6 mm ky −1 , respectively. Accounting for these low average erosion rates results in an average increase in sample exposure ages of~9%, which is a small modification to the exposure ages bounding the PBRs. We do not, therefore, modify PBR-bounding block removal exposure ages to account for this erosion for two reasons: first, because the effect on the exposure ages is minor and, second, because the average sea stack erosion rate is likely not spatially uniform at the scale of individual PBRs and therefore a systematic correction is not appropriate.
To determine the timing of block removal that established the current PBR geometry, that is, fragility age, we calculated the exposure ages of all samples collected from the seven PBRs at the Double Rock site. The exposure ages for the whole Double Rock site (n ¼ 31), and therefore timings of block removal, range from 17 ± 2 to 95 ± 9 ky (Figure 2e). A similar exposure age range for the whole Double Rock site is also observed at individual paleo-sea stack outcrop scale: 20 ± 5 to 95 ± 9 ky on DRE (n ¼ 10) and 17 ± 2 to 82 ± 8 ky on DRW (n ¼ 21). The individual PBRs show a similar age range around each ( Figure S4b): for example, 22 ± 2 to 39 ± 4 ky on DRW1 (Figure 2d). In summary, the distribution of surface exposure ages is similar at three different spatial scales of observation: at the whole Double Rock site (n ¼ 31, all 7 PBRs), each paleo-sea stack rock outcrops (n ¼ 10 and 21, 2 and 5 PBRs, on DRE and DRW, respectively), and each individual PBR (n ¼ 2 to 6 on each PBR). Based on these similar age distributions at different spatial scales, we interpret that all seven PBRs share a common evolution process in space and time, and we thus use the distribution of block removal exposure ages from all seven PBRs at the whole Double Rock site to interpret a single site-wide fragility age.
This distribution of exposure ages at the whole Double Rock site provides further support to our conceptual model of PBR evolution and fragility development by block removal events through time. The age at which the PBR's current geometry and therefore fragility, that is, fragility age, was achieved is recorded by the youngest exposure age bounding the PBR. The aggregated probability density function (PDF) considering all samples for all seven PBRs includes six samples generating the youngest peak ( Figure 2e); this youngest cluster of exposure ages gives an average fragility age of 21.0 ± 2.4 ky (90% confidence interval, 1.65σ; reduced χ 2 ¼ 0.79, P ¼ 0.56), which we interpret as the best estimate fragility age for all seven PBRs (section 3). In other words, all of the PBRs have existed in their current geometries for the last 21.0 ± 2.4 ky, with negligible change to their fragilities between that time and now. This fragility age also suggests that sufficient time has elapsed for the PBRs to have experienced strong ground motions from multiple moderate to major earthquakes on one or more of the local faults.

Validation of PSHA Hazard Curves
The fragility of each PBR and site-wide fragility age are then used in combination to place constraints on and reduce uncertainties in the seismic hazard curves from the 2015 DCPP PSHA model modified to use a non-ergodic GMM (Abrahamson et al., 2019;PG&E, 2015a) and computed for the Double Rock site (section 5.5.1). A probabilistic evaluation of the unexceeded ground motions interpreted for the PBRs at Double Rock is conducted using a method developed to test these site-specific hazard estimates based on PBR fragility functions Hanks et al., 2013). These fragility functions are the product of the site's PSHA vector hazard with the probability of toppling of each PBR to compute the rate of toppling for each PBR (Figure 3a and section 5.5.2). First, we assess the potential constraints that each of the seven Double Rock PBRs could place on the hazard curves. Specifically, we determine the scale factor that would need to be applied to the mean hazard curve in order for it to be consistent with the preservation of each PBR over the fragility age of 21 ky at a 5% probability of survival (β 0.05 ). The β 0.05 scaling factor of the mean hazard curve was calculated (section 5.5.3) to achieve the 5% target probability of survival of each PBR (Table S2).
In order to visualize the constraints placed on the fractile hazard curves by each PBR, the mean hazard curve scaled by the β 0.05 of each PBR is plotted along with the uncertainty fractile hazard curves. The scaled hazard curves are only shown at ground-motion amplitudes contributing to 25-75% of the failure probability for each of the PBRs (Figures 3a and 3b), which reflects the ground-motion range best constrained by the PBRs. For example, for DRW1, the PGA values that lead to 25-75% of the probability of failure range from 0.34 to 0.64 g. The scaled hazard curves are, therefore, parallel to the portion of the mean hazard curve over this ground-motion range. The fractile hazard curves are deemed to be consistent with the survival of the PBRs in their untoppled states if they either intersect or lie below the β 0.05 -scaled mean hazard curve over the range of ground-motion amplitudes representing the 25-75% failure contribution. Conversely, the hazard curves are regarded as being inconsistent with the PBR observational data if the fractile hazard curves estimate ground motions greater than the β 0.05 -scaled 25-75% failure ground motions of a PBR. For example, if the scaled hazard curve falls near the 60th fractile hazard curve, then 40% of the weighted hazard curves from the PSHA logic trees are regarded as being inconsistent (less than 5% chance) with the PBR not having been toppled.
At our study site, all seven of the PBRs demonstrate the potential reduction in ground-motion uncertainties by identifying inconsistent fractile hazard curves. Logically, it is the most fragile PBR, which in this study is DRW1, that provides the greatest constraint on the PSHA model ( Figure 3b). The β 0.05 -scaled mean hazard curve for DRW1 falls near the seventh fractile hazard curve; therefore, assuming that the individual hazard curves do not cross (see section 5.5.3), 93% of the weighted hazard curves are inconsistent with DRW1 in its untoppled state. We, therefore, use the most fragile PBR (section 3), DRW1, to place constraints on and reject the inconsistent hazard curves for the Double Rock site.
Second, to implement these PBR constraints and reduce uncertainties in the PSHA model estimates, we use the median fragility function and best estimate fragility age to identify the individual PSHA hazard curves for which there is a greater than 95% probability of failure (less than 5% probability of survival) of DRW1 (Equations 5 and 6 and section 5.5.4). The logic-tree end branches with hazard curves that are inconsistent with the PBR preservation are then removed from the PSHA model; that is, their end-branch weights are set to zero. The weights of the remaining logic-tree end branches that are consistent with the PBR preservation are renormalized so that the sum of all weights on all logic-tree branches again sums to 1.0 (Benjamin, 1970), and a new suite of hazard curve fractiles are calculated from these renormalized and reweighted branches ( Figure 4b). The revised PSHA results, which use the constraints of our best estimate fragility (fragility function median of 0.45 g, Figure 3a) and fragility age (21.0 ka, Figure 2e) of DRW1, have a new, PBR-informed mean hazard curve with ground motion levels that are reduced from the original mean curve. At an annual rate of exceedance of 10 −4 yr −1 , which is relevant for the design of some critical infrastructure, the level of ground motion associated with the PBR-informed mean hazard curve is reduced by 27% from the original PSHA model, from 1.02 to 0.74 g, following consideration of the PBR constraints. In addition to the change in the mean hazard, the uncertainties in the PBR-informed PSHA model are also reduced, as is evident from the reduction in spread of output hazard curves (epistemic uncertainty), and therefore increased precision, as a result of the PBR constraints ( Figure 4a). The range in the PGA for the 5th-95th fractiles at an annual rate of exceedance of 10 −4 yr −1 is reduced by 49%, from 1.04 g (before PBR constraints: 5th fractile ¼ 0.47 g and 95th fractile ¼ 1.51 g) to 0.53 g (after PBR constraints: 5th fractile ¼ 0.44 g and 95th fractile ¼ 0.97 g) and, in other words, is reduced by from a factor of 3.2 (0.47 to 1.51 g) to a factor of 2.2 (0.44 to 0.97 g).

Discussion
The fragilities and fragility age determined for the PBRs at Double Rock provide previously unavailable empirical data to analytically validate and update the PSHA model hazard results, which enables us to significantly reduce the uncertainties in the hazard curves output from the model at the site by 49%. It is important to note that PBR constraints provide an upper limit on PSHA hazard estimates, which results in hazard curves from the higher fractiles being removed from the PSHA model. Therefore, the application of PBR constraints will always result in a reduction in the hazard estimate of the site. This discussion section will first address the validity of the assumptions made in our exposure age calculations and, second, the validity of our fragility age. We will then justify the selection of the PBR probability of survival used in our hazard model evaluation and the decision to use only the constraints from the most fragile rock to remove hazard curves inconsistent with the PBR constraints at the site. Finally, we will establish the utility of our method for the identification of logic-tree branches in the PSHA model that are inconsistent with the PBR empirical data.

Accuracy of Individual Exposure Ages
The accuracy of the cosmogenic isotope surface exposure ages (Gosse & Phillips, 2001) calculated by our methods is evaluated by not only determining the rate of erosion of the chert sea stack outcrops both before and after the formation of the PBRs, but also the sensitivity of our ages to the assumption of block and outcrop self-shielding prior to surface exposure by block removal. First, our exposure age calculations are based on our geomorphic model, which assumes high erosion rates in the surf and splash zone prior to sea stack abandonment and PBR formation. The extremely low 10 Be concentrations measured in the three samples collected from the modern wave-cut notch are consistent with rapid erosion of the sea stack when subjected to wave action. This rapid erosion is removing previously produced cosmogenic 10 Be, as well as making it nearly impossible to preserve any PBRs that form at shore level for long periods of time. The inability of PBRs to be preserved while the sea stack is in the wave-action zone supports our geomorphic model of PBRs only being preserved following sea-level fall and associated platform abandonment. Therefore, we argue that the age of abandonment of the marine platform is a robust maximum age of fragility for all the PBRs at the site.
Second, in order to test the assumption of negligible erosion of the PBRs and surrounding surfaces since their exposure, we quantified the erosion rates of the paleo-sea stacks, and therefore the PBRs preserved on them, following the sea-level fall and abandonment of the sea stacks. Our calculated low erosion rates (2.4-2.6 mm ky −1 ) support the geomorphic model of the stability of the landscape following the abandonment of the major alluvial fan complex. These low erosion rates also support our interpretation that only minor modification has occurred to the PBR geometries that would affect the fragility age or fragility estimates we make. We therefore argue that the effect of erosion is minimal in our surface exposure age calculations (~9%). Furthermore, most of the colluvium volume used in the erosion rate calculation (section 5.3) is composed of large blocks that have fallen from the outcrops, similar to DRB-2 (height ¼ 4.6 m, width ¼ 1.6 m, depth ¼ 0.9 m; Figure S5), and this erosion is accounted for in the bounding ages of the PBRs formed by the removal of these large blocks. Therefore, the millimeter-or centimeter-scale erosion and resulting erosion depths become even more negligible in modifying both the surface exposure ages and the geometry of the PBRs. Finally, by using the minimum age of fan surface abandonment in the erosion rate calculations, our results yield a maximum rate of erosion; therefore, we are likely to be overestimating, rather than underestimating, the effect of erosion.
Third, we tested the assumption of complete block self-shielding prior to surface exposure by block removal. We calculated the shielding factor for a theoretical sample located on a vertical outcrop surface behind and shielded by a block of site-average height (2.56 m) and site-average width (2.13 m) and varied the block depth based on the range of observed block depths (1.10-0.38 m) at the site (section 5.4). Using these input parameters and a site-average block depth of 0.68 m results in a calculated shielding factor of 0.073 ( Figure S3), which means that the cosmogenic nuclide production rate at this theoretical sample point on the vertical outcrop surface behind this block of site-averaged dimensions is only 7.3% of that for an unshielded horizontal surface. In other words, our surface exposure ages assume that each surface that we sampled was 100% shielded by such a block and outcrop prior to block removal; however, we estimate that, for a theoretical sample point behind an average block at our site, only 92.7% of the incoming cosmic ray flux is shielded by the rock and outcrop prior to block removal. This difference between the assumed and estimated shielding (i.e.,~7.3%) using the average block depth has a relatively minor impact on our exposure ages and is similar to our fragility age uncertainty itself. Fortunately, our depth of measured blocks is likely biased toward relatively small depths because we favored fragile blocks that inherently have high height/depth ratios. Therefore, it is likely that the true site-average block depth for the rock outcrops is higher than the average value of 0.68 m and that our samples were even more completely self-shielded prior to exposure by block removal. For example, when the block thickness is greater than 0.92 m, which is within our observed range of measured block depths, the shielding factor is less than 0.05, which is less than the uncertainty of our fragility age. In summary, we argue that our incomplete knowledge of the shielding geometry of our samples prior to exposure by block removal has a negligible effect on our fragility age or associated conclusions.

Age Population Distribution and Selection of Fragility Age
For our study site, we implemented a new geomorphic model of PBR formation in order to establish a fragility age of the paleo-sea stack PBRs. We performed exposure age dating on samples collected from multiple bounding surfaces of multiple PBRs at the Double Rock site; this sampling strategy allowed us to use the population of bounding surface exposure ages to resolve the rate and timing of erosional processes acting on the paleo-sea stacks and interpret the PBR's fragility age. Our method not only allows the reconstruction of the timing of block removal events, but also it sheds light on the cause of these events. If the timing of erosion by block removal at the site occurred primarily as a result of discrete climatic or tectonic events, we would expect that the population of block removal ages at the site would cluster in time. We observe, however, weak clustering of ages ( Figure 2e); in fact, the population of the block removal ages at the site is generally quasi-randomly distributed in time at three different spatial scales of observation (whole Double Rock site, each paleo-sea stack, and individual PBR). Thus, we interpret that the seven PBRs included in our study at the Double Rock site share a common evolution process by quasi-random stochastic block removal prior to~21 ka.
An important exception to this quasi-random exposure age distribution is a clustering of exposure ages at 21 ky and no further block removal over the past 17 ± 2 ky ( Figure 2e); the selection of a best estimate fragility age of the PBRs is based on our observation of a peak in the PDF of all sample exposure ages. This peak at~21 ky, produced by the exposure ages of six samples, and, therefore, by the removal of six blocks, could be attributed to the last seismic event of sufficient magnitude to shake down fragile blocks from the outcrops. However, the lack of block removal ages younger than~17 ka could alternatively be correlated to the transition to a warmer, drier climate at the end of the Last Glacial Maximum (LGM), which may have resulted in a reduced surface erosion rate and thus increased preservation of the outcrops. We cannot deduce the cause of this hiatus in block removal since 17 ± 2 ky; however, given the high resistance to weathering of the chert outcrops in any climate, we speculate that it is more likely that the hiatus in block removal since 17 ± 2 ky reflects a lack of earthquake shaking of sufficient magnitude to shake down fragile blocks from the outcrops. Regardless of the cause of this cluster of exposure ages and associated peak in the exposure age distribution, the change in the rate and timing of erosion at~21 ka is interpreted as a transition from a period of PBR formation by block removal to a period of stability and PBR preservation. We, therefore, select 21 ± 2.4 ka (see section 2.2) as our best estimate fragility age of the seven PBRs studied at the Double Rock site when validating the consistency of hazard curves with the PBR constraints.

Selection of PBR Constraints
Our new method of refining PSHA models benefits from inherent versatilities that can be adapted to future PBR studies. For example, when selecting the constraints to be placed on the hazard, a judgment call can be made as to the selection of the PBR survival probability and the number of PBRs to use in the overall evaluation. Previous studies generally did not investigate both the distribution of fragilities and fragility ages of PBRs preserved at the site (Brune, 1996;. Our approach of extensively surveying a site for PBRs and analyzing a site-wide population of high-quality PBRs provides more empirical data to support the unexceeded ground-motion constraints that PBRs provide at the site. At the Double Rock site, the PBR data have been improved by obtaining a set of fragilities for seven PBRs and using the exposure ages surrounding all seven PBRs to interpret the process of their formation and, in turn, a site-wide fragility age. Such a population of fragilities and exposure ages facilitates a more informed selection of the PBR data that provide the most robust constraint on the hazard. The first versatility in our approach is the selection of the PBR survival probability to use to validate the consistency or inconsistency of PSHA output curves with the survival of the PBRs (Equations 5 and 6). We chose to use a threshold of 5% probability of survival (95% probability of failure) to be consistent with traditional levels used for hypothesis testing Hanks et al., 2013). We believe this to be an appropriate probability of survival for our study, given that we use the constraints of only one PBR from our analyzed population. We suggest that our 5% probability of survival provides not only a reduced risk of inaccurately removing PSHA hazard curves compared to a less conservative probability of survival, but also provides the ability to place significant constraints on the PSHA hazard curves. In future studies, the selected probability of survival can be adapted as appropriate to site characteristics and the requirements of the PSHA analysis.

AGU Advances
The second versatility in our approach is the selection of the number of PBRs in the analyzed population to use as constraints on the PSHA model. Considering our full set of fragilities from the seven studied PBRs at the site, and because we interpret all PBRs at the site to have the same fragility age, we chose to use the constraints of the most fragile rock, DRW1. In principle, the use of two rocks, DRW1 and the second most fragile PBR, DRW3, in the analysis would reject the hazard curves that are inconsistent with the two rocks and would therefore reject these PSHA outputs based on multiple observations. In practice, however, if both PBRs are required to topple during the time covered by their fragility age, then the site's fragility used to place constraints on the hazard is defined only by the fragility of the least fragile rock of the two, that is, DRW3. Although constraining the hazard with DRW1 and DRW3 would be doing so using multiple observations, compared to constraining the hazard with DRW1 only, we have the same confidence in the geomorphic, geologic, and tectonic assessment of DRW1 as the rest of the PBRs in our analyzed population. Therefore, we use only DRW1 as the maximum ground-motion constraint over the fragility age of~21 ky. However, in future PBR studies, alternative decisions on the number of PBRs to include can be made within our methods and adapted to the distribution of fragilities and fragility ages at that site.
Finally, we conducted a sensitivity analysis to assess the impacts of the uncertainty in the fragility and fragility age of DRW1 on the hazard curve constraints (Figure 4b). The epistemic uncertainty on the fragility age is estimated using the 90% confidence interval (±1.65σ) of ±2.4 ka calculated on the population of samples included in the mean 21 ka fragility age (Figure 2e). The epistemic uncertainty on the PBR fragility, however, is more difficult to estimate. Recognizing this limitation, and to allow for uncertainties in the process and determination of the PBR fragility, for example, the geometry and properties of the PBR-pedestal contact (Veeraraghavan et al., 2017), we assumed that the uncertainty on the median fragility is log-normally distributed with a standard deviation of 0.3 natural log units. With this assumption, the 5th and 95th percentile median fragility for DRW1 correspond to 0.27 and 0.74 g, respectively. The resulting nine different fragility and fragility age combinations all produce constraints that lead to the new normalized mean hazard curve being at a reduced ground-motion level relative to the original mean hazard curve. It logically follows that the lower-bound confidence level fragility (0.27 g) and the upper-bound confidence level fragility age (23.4 ka) produce the most significant constraint on the hazard curves. Conversely, the upper-bound fragility (0.74 g) and the lower-bound fragility age (18.6 ka) produce the least significant constraint. However, even with these conservative assumptions, there is still a significant reduction in the mean hazard. For example, the PBR constraints in these two end-member scenarios result in a reduction of the ground-motion level associated with the mean hazard curve at 10 −4 yr −1 of between 8% and 51% (0.87 to 1.65 g). Using the rigorous method with which we determine a fragility age and its uncertainty, the fragility age for DRW1 is well constrained, and the remaining small epistemic uncertainty in the fragility ages has a small effect on the hazard curves ( Figure 4b). If previous, less rigorous, methods were used, then the uncertainty in the fragility age would be much greater and reduce the robustness of the PBR constraints on the hazard.
The selection of fragility model also clearly impacts the constraints on the hazard curves. We showed that fragility estimates based on deterministic dynamic toppling thresholds are not best suited for the evaluation of a PSHA hazard model. For example, for DRW1, our estimate of the deterministic dynamic toppling acceleration is 0.23 g, whereas our probabilistic fragility function for DRW1 estimates a median PGA for failure of 0.45 g with a 5th-95th percentile uncertainty range of 0.27-0.74 g. The deterministic dynamic toppling acceleration is near, but below, the fifth percentile from the probabilistic results, which will underestimate the probability of survival. Additionally, the deterministic dynamic toppling value of 0.23 g falls outside the ground-motion levels best constrained by PBR 25-75% probability of failure ground motions of 0.34-0.64 g. Our study provides striking evidence that conventional deterministic fragility methods can significantly underestimate the ground-motion levels required for PBR toppling. This mischaracterization of the PBR fragility, in turn, could result in inaccurate evaluation of the PSHA model and the possibility of erroneous rejection of hazard curves. Therefore, probabilistic fragility methods should be used in PBR studies when evaluating PSHA models.

PSHA Model Weight Renormalization
To improve the PSHA model, we renormalized the model to contain only the hazard curves consistent with the ground-motion constraints provided with the PBR, DRW1. We achieved this by first setting the branch end weights of these inconsistent branches to zero. Then, the PBR-consistent logic-tree branch end weights 10.1029/2020AV000182

AGU Advances
were renormalized, ignoring the original node of the logic tree, to make their total sum to unity again. This renormalization is needed because the sum of all weights over all logic-tree branches must equal 1.0, which was the case before any branches were rejected and must remain true after any branches are removed. Using this renormalization method, we show the utility of characterizing the fragility and fragility age of PBRs to eliminate hazard curves that are inconsistent with the PBRs and thereby reduce uncertainties in hazard estimates, but it is beyond the scope of this paper to present the most appropriate method to renormalize and revise the branch weights for individual logic-tree nodes in the PSHA model.
Follow-up research could test and compare alternative probability methods to perform the reweighting of the logic-tree branches. Such research could statistically analyze the branches of each node in the logic tree that are associated with the rejected end-branch hazard curves and, in theory, identify which alternative input parameters and models (e.g., fault slip rates, magnitude-frequency distributions, and GMMs) should be targeted for re-evaluation and improvement in future studies. Automatic methods, such as Bayesian updating, could be trialed to revise the branch weights on such individual nodes. The advantage of such methods would be that the decisions regarding the choice of cutoff probability are avoided. The application of Bayesian updating to give new logic-tree branch weights is straightforward from a mechanistic perspective; that is, the initial distribution of hazard curves represents the prior distribution, and the likelihoods of the PBR survival can be combined with this prior to obtain the posterior distribution; however, there are conceptual nuances that warrant further consideration when using such purely statistical approaches. For example, logic-tree branch end weights associated with each hazard curves result from the combination of branch weights on individual nodes across the logic tree. A branch end weight associated with a hazard curve could have its weights modified via Bayesian updating, but this modification to the branch end weight would also effectively modify each of the branch weights on the individual nodes associated with that hazard curve. Such a modification of branch weights on individual nodes could result in the same branch of a logic-tree node having different effective weights depending on how it combines with downstream branches. On one hand, these different effective weights are not a problem when each path through the logic tree and the associated weights and branches are viewed collectively as a set of decisions or parameters associated with the hazard calculations. However, on the other hand, it could be a problem when branch sets at logic-tree nodes are developed largely independently from the rest of the tree, which is commonly done within hazard analyses. Reconciling these issues is an ongoing area of investigation, and trialing automatic methods was deemed beyond the scope of the present study. We recommend that, based on the judgment of experts performing the PSHA evaluation, the hazard curves that are either rejected or whose weights are heavily down-weighted following Bayesian updating are carefully inspected to understand which logic-tree branch combinations led to those hazard curves. If particular logic-tree branches are responsible for the sets of hazard curves that are rejected or down-weighted, then one would have a strong argument to reassess the validity of those particular branches. Alternatively, rather than just change the weights on the existing branches, the evaluation of the individual nodes may conclude that the current models are inadequate and new models need to be developed and added to the logic tree.

Conclusions
Our new method provides previously unavailable probabilistic constraints to quantitatively evaluate a site's PSHA model and reduce the uncertainties in the hazard estimate by rigorously assessing the fragilities and fragility age of a population of PBRs. The PSHA curves that are inconsistent with the PBR observations at the threshold probability level can be removed and the remaining end branches renormalized to provide revised, improved PSHA model results. At ground-motion levels with an annual rate of exceedance of 10 −4 yr −1 at the Double Rock site, we not only reduce the mean ground-motion estimate by 27% but also reduce the range of the 5th-95th hazard fractiles by 49%. The Double Rock site in Central California is an ideal example of the application of this new method because PBRs have formed proximal to important seismic sources, where there is sparse ground-motion data to constrain the non-ergodic part of the GMM, and at a site near critical infrastructure, the Diablo Canyon nuclear Power Plant (DCPP), which has recently been subjected to a high-level (SHACC level 3) PSHA study that followed what is considered the current best international practice (SSHAC, 2012). Specifically, the Double Rock site is located at a similar distance from the Hosgri fault zone and on similar site characteristics as DCPP. Therefore, our work tests the ground motions resulting from the same seismic sources that control the hazard at the DCPP site.

10.1029/2020AV000182
Our study broadens the utility of PBR studies globally to address fundamental and societally relevant issues in earthquake science. First, we show that the PBRs provide constraints on the PSHA hazard model over a timescale of 21.0 ± 2.4 ky. The preservation of the PBRs over this timescale provides uniquely valuable constraints on the rates of rare, large-magnitude earthquakes and their ground motions, which have large uncertainties because few such events have been recorded by modern instrumentation. The vast majority of previous PBR studies (Brune, 1996; have focused on granitic terrains in southern California and commonly assume a default fragility age of 10 ka if PBR-specific age constraints were not obtained (Grant Ludwig et al., 2015). Our study shows that PBRs can be preserved in the landscape for twice as long as conventionally assumed and, therefore, potentially provide ground-motion constraints over significantly longer timescales. Furthermore, this study is the first analysis of PBRs in a coastal geomorphic setting and where the PBRs formed by the process of fracture-bound block removal on paleo-sea stacks located on tectonically uplifted marine terraces. Uplifted marine terraces are common geomorphic features along tectonically active coastlines throughout the world, many of which are likely to contain paleo-sea stacks with optimally oriented rock fabric sets to form PBRs (e.g., along the coast of the Pacific Northwest, USA, and the Cascadia subduction zone). We have developed a new geomorphic model for how such PBRs can form and a sampling strategy to determine their fragility age. Thus, our techniques and this method have the potential to be used broadly to test and update site-specific hazard curves for coastal areas including regions in which the controlling seismic sources are offshore faults that are inherently more difficult to constrain. Beyond coastal regions, our method of PSHA result revision and reduction in uncertainties is generally applicable to site-specific hazard studies globally to mitigate the ubiquitous lack of near-field, large-magnitude ground-motion recordings at the site of interest. The empirical geologic data for validating and refining hazard estimates delivered by our method are essential to reducing the very large uncertainties in site-specific hazard estimates so that the results can be reliably used in decision making regarding adequate seismic safety in the numerous earthquake-prone areas where critical infrastructure is located near human populations.

3-D Model Construction
For the seven PBRs at the Double Rock site, a set of photographs that captured all available view angles was collected using a calibrated hand-held camera. The commercial photogrammetric software PhotoModeler was used to co-locate common points among photos to orientate the photos relative to each other. The iterative solver within PhotoModeler generates 3-D points within a~1% accuracy of their true locations (Anooshehpoor et al., 2013). The geometry of the PBRs was manually constructed, within a defined coordinate system, as points defining planar surfaces on the PBR. The complex geometry of the PBRs caused auto surfacing tools within PhotoModeler to produce inaccurate results; therefore, the points were triangulated manually among each other to produce a surfaced 3-D model. These PhotoModeler triangulated models have been shown to have volumes within a few percent of the objects' true volumes (Anooshehpoor et al., 2013).

Geometric Parameters
The 3-D models were then imported into MATLAB to calculate the center of mass, and rocking points were selected along the PBR-pedestal contact in both fragile directions. The Double Rock PBRs, like conventional granitic PBRs studied in southern California (Brune, 1996;, are free to oscillate about the rocking points in both fragile directions. However, unlike conventional PBRs, the outcrop surface on one side of each of the PBRs dictates that the PBRs are only free to topple in a single direction, that is, away from the outcrop, and are, therefore, more similar to some of the PBRs studied at Yucca Mountain ). An additional similarity to some of the PBRs studied at Yucca Mountain is the unconventional multi-block geometry of two of the Double Rock PBRs (DRW4 and DRW5); the geometry of each of these multi-block PBRs was simplified and modeled as a single coherent block.
Toppling of a PBR is controlled by two geometric parameters (Purvance, Anooshehpoor, et al., 2008): a radius (R) that connects the center of mass of the PBR and the rocking point with the angle (α) between the radius and vertical (Figures 2a and S2a). The critical rocking point for toppling is defined by the point with the smallest value of α (α min ). Three iterations of critical rocking point searches were conducted for each PBR, with each subsequent iteration occurring between the two points that previously produced the lowest α values. The α min angles in both fragile directions were compared to verify the symmetric geometry and associated rocking response of each PBR. For a symmetric PBR, α 1 ¼ α 2 , whereas for an asymmetric PBR, α 1 < α 2 . A difference in the two α values of less than a 5% was used to classify the PBRs as symmetric and indicated an approximately horizontal pedestal.
At the Double Rock site, the two most fragile PBRs (DRW1 and DRW3) have symmetric geometric factors (Table S1); that is, α 1 is within 5% of α 2 , and, therefore, it is straightforward to use the appropriate symmetric rocking equations of the fragility model to calculate the fragility of these symmetric PBRs (Purvance, Anooshehpoor, et al., 2008). The other five analyzed PBRs at the Double Rock site, however, are less straightforward because they have asymmetric geometric factors (Table S1); that is, α 1 and α 2 differ by more than 5%. Although the fragility model also includes equations to describe the rocking response of PBRs with asymmetric geometric factors, these rocking equations for asymmetrical PBRs assume that α 1 < α 2 and that the PBR will topple in the direction of α 1 (Purvance, Anooshehpoor, et al., 2008); however, this assumption is violated for the asymmetric PBRs at the Double Rock site because their geometry is such that they cannot topple in the direction of α 1 because α 1 is toward the outcrop surface. In other words, the asymmetric PBRs at the site can only topple in the α 2 direction, that is, away from the outcrop, and therefore, the asymmetric rocking equations of the fragility model are not appropriate for these asymmetric PBRs.
Considering the lack of rocking equations appropriate to the asymmetric PBRs at the site, we instead made a conservative best estimate of the fragility for these PBRs by using the symmetric rocking equations with the α and R calculated in the only direction that the PBR can topple, that is, using α 2 and associated R 2 in the direction away from the outcrop. There are two important reasons why this method provides a conservative estimate of the fragility for these asymmetric PBRs at the site. First, asymmetric PBRs will always be less fragile than symmetric PBRs with the same minimum α 1 because rocking motion is dampened in the larger α 2 direction of asymmetric PBRs relative to the rocking motion of symmetric PBRs. Second, these Double Rock asymmetric PBRs always have a greater α and, therefore, are less fragile in the direction that they can topple away from the outcrop (Table S1).
In using this method for the asymmetric PBRs at the site, we acknowledge that we may be mischaracterizing the true rocking response of the PBRs with asymmetric geometric factors; however, any mischaracterization is likely to underestimate, rather than overestimate, the fragility of the asymmetric PBRs at the site. Furthermore, it is important to note that neither our validation of the PSHA hazard curves results nor any of our conclusions derived from them are sensitive to such a mischaracterization because only the most fragile, symmetrical PBR with straightforward geometric factors is used to inform our subsequent hazard model evaluation. Specifically, DRW1, which is the PBR used to modify the PSHA model, has a simple symmetric geometry. We are, therefore, confident that the rocking equations used for the most crucial PBR, DRW1, are entirely appropriate to estimate its fragility. That said, future work, such as additional shake table experiments or finite-discrete element modeling (Veeraraghavan et al., 2017), could improve the accuracy of the fragility analyses of asymmetric PBRs such as those at the Double Rock site, for example, that can only topple in one direction toward α 2 and with multi-block PBR geometries.
The slope of the pedestal was shown to be a primary control in the direction of PBR toppling rather than the directionality of three-component acceleration histories (Veeraraghavan et al., 2017). Experiments controlling the direction of initiation of PBR rocking showed that the rocking could be initiated in any direction, but the PBR then finds it easiest to topple in the vicinity of the minimum α due to the absence of strong directionality in the ground-motion recordings (Veeraraghavan et al., 2017). The accurate assessment of the contact geometry and conditions between the PBR and its pedestal are, therefore, required for accurate fragility analyses. However, these are very difficult to investigate without removing the PBR and inspecting the geometry of the pedestal beneath the PBR. Such an invasive investigation would lead to permanent damage to the PBR. Alternatively, conducting field tests to determine the quasi-static toppling force of a precarious rock and frictional/cohesive forces at the PBR-pedestal contact (e.g., Anooshehpoor et al., 2004) was not feasible and beyond the scope of our study, primarily due to the mass of each of the PBRs (e.g., DRW1 ¼~10 tons) and their location high up on the outcrops. We, therefore, do not have information on the conditions of these contacts between each PBR and its pedestal.

Fragility
The tangent of the angle α min in the direction that the PBR can topple, that is, away from the outcrop surface, multiplied by the acceleration of gravity (g) approximates the minimum quasi-static horizontal acceleration threshold required to cause the PBR to topple (Weichert, 1994). A dynamic toppling acceleration threshold accounts for the nonlinear rocking response of a PBR and the dependence on the direction of maximum ground motion, which empirical data have shown to be, on average, 30% greater than the quasi-static acceleration (Anooshehpoor et al., 2004). This dynamic scale factor is applied to the PBRs to calculate a deterministic threshold ground acceleration that would cause toppling, given by Equation 1.

Deterministic Dynamic Toppling Acceleration
Empirical relationships regressed from shake table experiments determine the probability of toppling of a PBR based on its measured α min angle and length R value, for a PGA and PGV/PGA pair (Purvance, Anooshehpoor, et al., 2008). A PBR's fragility matrix is calculated by solving the empirical relationship (Equation A1 with corrected final coefficient in Purvance, Anooshehpoor, et al., 2008) for the probability of toppling given pairs of PGA and PGV/PGA values (Figures 2b and S2b), given by Equation 2. (2)

Sample Collection
Samples were collected from rock surfaces surrounding all accessible sides of the seven PBRs to constrain the timing of removal of the surrounding blocks that created the PBRs' current shape. Continuous outcrop surfaces that suggested the past location of large blocks were selected for sampling; these surfaces were highly varnished, which suggests that they had not been recently superficially modified by erosion after block removal. Samples were preferentially collected from vertical outcrop surfaces that we interpreted were originally located behind the base and center of these large blocks in order to maximize the amount of shielding of the sample prior to its exposure by block removal. Additionally, a sample was collected from a detached block that fell onto the major fan surface (DRB-2, Figure S5), and three samples from a modern wave-cut notch in the intertidal zone ( Figure S1) were also processed by the following methods. Approximately 1 kg of chert bedrock material was collected per sample, using hammer and chisel, and its thickness recorded.

Sample Preparation
The majority of sample preparation was conducted in the CosmIC Laboratory at Imperial College London. Samples for the PBR DRE1 were prepared at the Scottish Universities Environmental Research Centre (SUERC). Samples were crushed to 250-500 μm grains and purified to SiO 2 mineral separates using magnetic separation followed by acid etching to remove any atmospherically derived 10 Be adhered to the grain surfaces. Beryllium was isolated by anion and cation exchange and prepared into targets for AMS analysis (Corbett et al., 2016).

10
Be samples were primarily analyzed at the Australian Nuclear Science and Technology Organisation (ANSTO) using the 6MV SIRIUS accelerator mass spectrometer (Wilcken et al., 2017) during analysis runs in April and June 2018. Samples associated with the PBR DRE1 were analyzed at SUERC in May 2014 (Xu et al., 2015). 10 Be/ 9 Be data from ANSTO were normalized to the primary standard KN-5-2 with a nominal value of 8.558 × 10 −12 (Nishiizumi et al., 2007), whereas data from SUERC were normalized to NIST standard with an assumed value of 2.79 × 10 −11 (Nishiizumi et al., 2007); in both AMS laboratories, two secondary 10.1029/2020AV000182 AGU Advances standards were run as unknowns to confirm the linearity, accuracy, and precision of the measurements. The measured 10 Be/ 9 Be ratios were reduced to the number of total 10 Be atoms in the sample using the mass of low-background beryllium carrier added to each sample (UVM 17-1 beryl carrier with~3 × 10 -16 10 Be/ 9 Be). A process blank, of beryllium carrier only, was processed and analyzed with each batch of samples.
The average calculated number of total 10 Be atoms for all measured batch blanks in each accelerator run (Table S2) was subtracted from every sample's calculated number of 10 Be atoms, and each sample's 1σ AMS analytical error and the standard deviation of process blanks in that run were propagated in quadrature.
The number of 10 Be atoms (N 10 ) and uncertainties was reduced to 10 Be concentrations and uncertainties using each sample's mass of quartz dissolved, following Equation 3.
where M q is the mass of the sample quartz, R 10/9 is the AMS measured ratio of 10 Be to 9 Be, M c is the mass of 9 Be added as carrier, N a is Avogadro's number, A Be is the molar weight of Be, and n 10. B is the number of 10 Be atoms in the process blank. All the necessary data reduction inputs to be used in Equation 3 are provided in Tables S4 and S5. A replicate sample (DRW1-D) was processed for both the April and June 2018 ANSTO accelerator runs; the measured 10 Be concentrations for the two samples were within 1σ uncertainty.

Shielding Calculations
The complex outcrop geometry obstructing the cosmic ray flux to each sample point, which modulates the 10 Be production rate in each sample, was individually calculated from 3-D models created by photogrammetry of the PBR and its adjacent outcrop. Sensitivity analyses were conducted, and the outcrop models extended a minimum radius of 2 m from the sample locations to ensure the accuracy of the calculated shielding factor. The code (Balco, 2014) uses Monte Carlo integration to model the relative attenuation of the cosmic ray flux by rock and air to each sample point, in order to calculate a dimensionless shielding factor relative to an unshielded sample. The Monte Carlo integrations were run for 1,000 iterations, with convergence occurring after~300 iterations. The calculated shielding factor to each sample is included in Table S6.

Exposure Age Calculations
The online CRONUS v2.3 surface exposure age calculator (Balco et al., 2008) was used to interpret each sample's exposure age and 1σ uncertainty (Figures 2c and S4) from its location, elevation, and thickness, as well as the calculated 10 Be concentration and shielding factor. Ages are calculated using a constant production rate model and "St" scaling scheme for spallation (Lal, 1991;Stone, 2000), with a 10 Be production rate of 4.132 ± 0.218 atoms g −1 yr −1 , based on the "primary" calibration data set of Borchers et al. (2016). Muonogenic 10 Be production is modeled after Heisinger et al. (2002aHeisinger et al. ( , 2002b. Additionally, a rock density of 2.7 g cm −3 , a standard atmosphere, and no erosion were assumed. The full CRONUS surface exposure age calculator inputs are provided in Table S6.

Erosion
An estimate of the long-term average erosion rate of the chert outcrops was made by calculating the volume of colluvium that has accumulated on the major fan surface since the fan's abandonment. To estimate the rate of erosion, we first need to estimate the age of the fan. This process is described further below.

Dating of Fan Surface Abandonment
The maximum age of the Qf 1 fan surface is the abandonment of the Marine Isotope Stage (MIS) 5a terrace, at 85 ka, assuming an immediate major depositional event following the sea-level fall. The minimum age of the fan surface has been constrained by the dating of a block (DRB-2, Figure S5), which fell from the DRE outcrop onto the abandoned fan surface, and thus postdates fan abandonment. A sample for cosmogenic exposure dating was collected from the rock surface that was originally maximally shielded by the depth of the block (~0.9 m) prior to it falling from the outcrop, that is, the reconstructed basal center of the block originally adjacent to the outcrop surface. Thus, the entire 10 Be inventory measured in the sample accumulated in the current fallen geometry ( Figure S5). Topographic shielding azimuths and elevations to the DRB-2 sample point were measured in the field using a compass and clinometer. The shielding factor was determined using the CRONUS online topographic shielding calculator (v.1). The online CRONUS v2.3 surface exposure age calculator was used to interpret the sample exposure age and 1σ uncertainty of the sample and thus the timing of deposition onto the abandoned major alluvial fan surface (Tables S5 and S6). This sample provides a minimum age for the major fan surface of 49 ± 4 ka, as the period of fan stability prior to the falling and deposition of DRB-2 is unknown.
The characteristics of the soil profiles developed at two locations at the Double Rock site were described and also used to estimate the length of time the alluvial fan deposits have been exposed to soil-forming processes at the ground surface (Figure 1b and Text S1; Dolan et al., 1997;Muhs, 1982;Rockwell et al., 1984Rockwell et al., , 1987Rockwell et al., , 1989Rockwell et al., , 1990Rockwell et al., , 1994Soil Survey Staff, 1999, 2014. The soil profiles were found to be comprised of a surface soil and several buried soils. The Buried Soil 1 in Profile 1 (DRP1) is thought to correlate with the Surface Soil in Profile 2 (DRP2). The Surface Soil in Profile 1 is interpreted to have developed on minor fan deposits that overlie the major alluvial fan deposits; these do not spatially extend to the location of Profile 2. These two soil profiles (including the surface soil and buried soils) have minimum estimates of the length of time of soil development of 56 and 49 ky, respectively. We use the estimated soil development ages and the 10 Be age calculated for DRB-2 (49 ± 4 ka) to infer a conservative (i.e., minimum) major fan surface age of 45 ky (mean age minus 1σ).

Colluvium Volume Calculation Method
A~1 km 2 terrestrial LiDAR survey of the site was conducted, and a 0.5 m gridded digital elevation model (DEM) was generated in ESRI ArcMap. The gradient of the unmodified major alluvial fan surface was computed, using natural neighbor interpolation, to extend beneath the colluvium and chert outcrops (Figure 1b). The colluvium extent around the southern faces of DRE and DRW was accurately mapped in ArcMap, and the volume difference between the interpolated fan surface and LiDAR DEM in this area was calculated.
The surface areas of the DRE and DRW chert outcrops, which are the source of colluvium, were determined using the terrestrial LiDAR survey point cloud and the software MeshLab. The point cloud was simplified and surfaced using a screened Poisson surface reconstruction. The model was clipped to the southern faces of the chert outcrops only, and the surface areas for each clipped outcrop were calculated. The colluvium volume (DRE ¼ 67 m 3 , DRW ¼ 153 m 3 ) was combined with the outcrop surface areas (DRE ¼ 610 m 3 , DRW ¼ 1,300 m 3 ) and the minimum fan age (45 ka) to calculate a maximum erosion rate (DRE ¼ 2.4 mm ky −1 , DRW ¼ 2.6 mm ky −1 ).

Block Self-Shielding Analysis
The average height, width, and depth of the eight blocks (seven PBRs and one fallen block) used in our study were measured on the PhotoModeler 3-D models constructed for each block. A theoretic block was then constructed in PhotoModeler with the average height (2.56 m) and width (2.13 m) of the blocks at the site ( Figure S3a), and a vertical outcrop of 10 m (height) by 12.13 m (width) by 5 m (depth) to ensure complete 1π shielding by the outcrop surface constructed behind the block ( Figure S3b). The depth of the theoretical block was varied across the range of measured block depths at the site. For each block depth, the shielding factor was calculated (Balco, 2014) for a theoretical sample located at the center of the bottom edge of the theoretical block where it met the outcrop surface ( Figure S3c). This theoretical sample is representative of the location of our surface exposure age samples, which we preferentially collected from vertical surfaces behind the base and center of where large blocks are believed to have been located before their removal exposed the sampled surface. The shielding factors for each block depth were plotted to show the decrease in shielding factor with increasing block thickness ( Figure S3d).

Validation of Seismic Hazard Estimates 5.5.1. Seismic Hazard Model
Typical PSHA utilizes GMMs that rely on the ergodic assumption and use average source, path, and site effects derived from global databases (Anderson & Brune, 1999). Our PSHA model for the Double Rock site, however, differs from a typical PSHA model because it uses a suite of alternative non-ergodic GMMs (Abrahamson et al., 2019) that reflect the systematic source, path, and site effects for a given earthquake scenario and a given site location. Ergodic GMMs model these systematic effects as randomness that can occur at any site from any scenario. As a result, ergodic GMMs have larger aleatory variability than is appropriate for a given site/source pair and will not represent the site/source-specific median ground motion (Anderson & Brune, 1999). It does not make sense to test hazard curves computed using ergodic GMMs with site-specific constraints on the hazard from the PBRs.
At the Double Rock site, there is sparse ground-motion data to constrain the site-and source-specific effects in the non-ergodic GMM, which leads to large epistemic uncertainties in the median GMM that are reflected in the broad spread of fractile hazard curves (Abrahamson et al., 2019). To reduce these large uncertainties, we need either ground-motion data from earthquakes near the Double Rock site or other constraints on the ground shaking at the site. While the PBRs do not provide a direct constraint on the GMM at the Double Rock site, the PBRs do provide constraints on the hazard curves, which indirectly provide constraints on the GMMs because the uncertainty in the GMMs tends to dominate the uncertainty in the hazard curves for long return periods. We, therefore, aim to reduce epistemic uncertainty in the hazard and the underlying non-ergodic GMMs by evaluating the constraints that the PBRs provide on the hazard curves.
The PSHA model we evaluate for the Double Rock site was computed using the 2015 DCPP PSHA model modified to use a non-ergodic GMM (Abrahamson et al., 2019;PG&E, 2015a) for site conditions of V s30 ¼ 1,000 m s −1 because of the hard-rock outcrops on which the PBRs are located and includes three alterative models for uplift based on reverse faults under the Irish Hills. Monte Carlo sampling of the logic-tree branches is used to generate a subset of 20,000 equally weighted hazard curves and provides a mean hazard curve, as well as a suite of fractile hazard curves. This PSHA model, on one hand, computes an estimate of scalar hazard for PGA ground-motion levels. The PBR fragility model (Purvance, Anooshehpoor, et al., 2008), on the other hand, computes a probability of toppling as a vector of both PGA and PGV/PGA ground-motion levels.
In order to evaluate the PSHA model using the PBR fragility model, we must first make the scalar and vector outputs of these models compatible with each other. In our case, we convert the PSHA model scalar hazard estimates from PGA ground-motion levels to vector hazard estimates, which describe the joint rate of exceeding ground-motion levels of both PGA and PGV/PGA derived from input time histories applicable to the site with a range of PGA-PGV/PGA combinations. To compute the vector hazard, the PGA hazard curve is transformed into a vector hazard using a conditional GMM for PGV that is based on the PGA (Equation 12 in Abrahamson & Bhasin, in review).

PBR Fragility Functions
The fragility model gives the probability of toppling the PBR for a given pair of PGA and PGV/PGA values. Multiplying this probability by the rate of occurrence of the PGA-PGV/PGA pair from the vector hazard, derived from input time histories with a range of PGA-PGV/PGA combinations, gives the rate of toppling of the PBR for this PGA-PGV/PGA ground-motion pair. The sum of these rates of toppling over all PGA and PGV/PGA pairs provides the annual probability of failure (P annual failure ) of the PBR. The vector fragility can be collapsed to a scalar fragility that depends only on PGA by summing the rates of toppling over the suite of PGV/PGA values for a given PGA value. This approach accounts for the PGV/PGA of the time histories that are applicable to the site. The cumulative contribution of all PGA ground-motion levels to the annual probability of failure (i.e., the probability of failure) is calculated, which we define as the fragility function of the PBR (Figure 3a) (cf. Baker et al., 2013;Hanks et al., 2013). The different combinations of PGV/PGA, as well as their relative likelihoods, that arise for different levels of PGA are directly accounted for when reducing the vector hazard space to scalar hazard space. The slopes of the fragility curves in Figure 3a reflect the variability of the seismic capacity of the PBRs, and the influence of PGV is reflected in that variability. A fragility surface, rather than a fragility curve (Figure 3a), would show steeper slopes in these fragility curves for a given value of PGV. From this fragility function, we can calculate the median 50th percentile probability of failure. We use this median ground motion and the slope of the fragility function as primary metrics used to describe the PBR fragility.

β Ps -Scaled PBR Constraints
To better visualize and quantify the potential constraints provided by the fragility function and fragility age of each of the seven PBRs, the scaled mean hazard curve that leads to a 5% probability of survival for each PBR is compared to the PSHA fractile hazard curves (Figure 3b). The β Ps -scaled mean hazard provides an estimate of the upper bound of the PSHA output hazard curves that are consistent with the existence of the PBR. Any PSHA output hazard curves that fall above the β Ps -scaled mean hazard curve would be rejected. For example, if the scaled mean hazard curve lies near the 80th fractile, that would indicate that the upper 20% of the logic-tree branch weights would be rejected. However, the β Ps -scaled mean hazard curve does not provide a direct estimate of the PBR-informed mean hazard curve as a result of inconsistent hazard curve rejection. The PBR-informed mean hazard curve must be recomputed after the inconsistent hazard curves are removed and the weights of the remaining consistent hazard curves are renormalized. It is worth noting that the renormalized hazard curves, as a result of removing even a single inconsistent hazard curve, will always have a mean hazard curve that is reduced from the original level.
The scale factor, β Ps , is calculated using where P annual failure is the probability of failure using the full set of logic-tree branches without considering the PBR constraints and T is the fragility age of the PBR. The value of P s is selected based on the confidence of the analyses used in defining T and P annual failure , as well as desired conservatism: Here, we use a P s value of 0.05 to calculate a scale factor β 0.05 Hanks et al., 2013). The value of β Ps describes the potential for the PBR to reject the upper fractiles of the hazard. To visualize the range of ground-motion levels for which the hazard is best constrained by the PBR, the PGA range that contributes between 25% and 75% of the total probability of failure from the PBR fragility function is shown ( Figure 3a); this 25-75% range of PGA is chosen to avoid constraining the hazard over ground-motion levels that are not well constrained by the PBR. To plot these PBR constraints in hazard space, the annual frequency of exceedance from the mean hazard curve associated with this 25-75% PGA range is scaled by β Ps to be plotted with the fractile hazard curves (Figure 3b).
In assessing which hazard curves would be rejected by the β Ps -scaled mean hazard curve, we implicitly assume that there is a clear ordering of individual hazard curves, that is, that the individual curves do not cross. We believe this assumption to be appropriate for the purpose of illustrating the potential reduction in hazard provided by each of our PBRs because the β Ps -scaled range shows the ground-motion level that contributes to the central half of the probability of failure, that is, the 25-75% range of PGA (Figure 3b). While hazard curves can cross over in this central range of the ground motions that contribute to failure, the main difference in the hazard curves will be a scale factor, which is represented by the range on the fractiles. As this range on the fractiles covers half the contribution to failure, this ground-motion level gives a good indication of the hazard curves that will be rejected. Additionally, the more accurate probability of failure is computed by our final analysis of each of the individual hazard curves (Figure 4).

Hazard Curve Rejection
The β Ps is useful for identifying the critical PBRs in a population and quantifying their potential constraints; however, it is not the best tool for identifying inconsistent hazard curves, as β Ps uses the shape of the mean hazard curve instead of the shape of each of the potential alternative curves. We improve upon previous methods Hanks et al., 2013) and separately assess each alternative hazard curve output by the PSHA model. To determine which of the hazard curves are inconsistent with the PBR observations, the entire PGA fragility function of a PBR (Figure 3a) can, in turn, be combined with each alternative PGA hazard curve output from the PSHA model. First, for each output PGA hazard curve, the rate of occurrence of PGA is obtained as the derivative of the hazard curve annual frequency of exceedance. Second, the rate of failure for the PBR, γ Failure , is obtained by combining the probability of failure for a given level of PGA, P(Failure| PGA ¼ a), with the rate of occurrence of that level of PGA, and then summing these combinations for all levels of PGA.
Using the approach described above and in Equation 5, the rate of failure (γ Failure ) of the PBR is used to calculate its probability of survival, P Survival in T , (i.e., no failure) over the fragility age, T, using Equation 6.
The tested hazard curve is rejected if the probability of survival over the fragility age exceeds a defined percentage: For this study, a threshold of 5% is used. The basis for selecting this threshold probability of survival follows statistical hypothesis testing convention Hanks et al., 2013). The hazard curves that are inconsistent with the PBR constraints are weighted to zero in the PSHA model, and the remaining logic-tree branch weights are renormalized to sum to unity (Benjamin, 1970). A new suite of hazard curves is then calculated from the subset of PBR-consistent hazard curves, from which new fractile hazard curves can be plotted (Figure 4a).
The 90% confidence intervals (±1.65σ) in both the fragility and fragility age of DRW1 are considered in our validation and rejection of the hazard curves. The fragility function median ground motion of DRW1, 0.45 g, which we estimate to have a log-normally distributed uncertainty, with a log-standard deviation of 0.3, has been used to investigate three levels of fragility. This log-normally distributed uncertainty is used instead of an epistemic uncertainty on the PBR fragility because an epistemic uncertainty cannot be calculated by our fragility method. The 5th and 95th percentile median fragility corresponds to 0.27 and 0.74 g. Three fragility ages (T) have also been considered: the best estimate 21.0 ka fragility age, and each limit of the 90% confidence interval (±1.65σ) on this age of ±2.4 ka. This allows for us to consider nine different potential impacts on the hazard curves using the PBR constraints (Figure 4b).