Bayesian analysis of ground motion models using chimney fragility curves: 2021, 5.9-Mw Woods Point intraplate earthquake, Victoria, Australia Open Access

Ground-motion models (GMMs) are a key component of probabilistic seismic hazard analysis (PSHA) (Bommer et al., 2010). In stable continental regions (SCRs) where moderate-to-large earthquakes are infrequent and the spatial density of seismic networks may be low, instrumental strong ground motion data can be rare, particularly in near-source areas (e.g. 0–50 km epicentral distances). This can compromise objective analysis of GMM parameters and performance in forecasting ground shaking intensities and distributions (Leonard et al., 2014).

In Australia, development of the National Seismic Hazard Assessment (NSHA) is led by Geoscience Australia (GA). The NSHA 2018 (Allen et al., 2018) used an expert elicitation process to select and weight GMMs used to develop the seismic hazard model (Griffin et al., 2018, 2020). An important aspect of selecting and ranking GMMs for use in hazard assessments is the comparison of instrumentally recorded ground motions with theoretical outputs from the GMMs (e.g. Ghasemi and Allen, 2018). The development of Australian-specific GMMs is a challenging task given the sparse recording networks and limited strong-ground-motion recordings. Three GMMs that have been specifically developed for the region are commonly used in Australian seismic hazard assessments: A12 (Allen, 2012), SEA09YC (Somerville et al., 2009—Yilgarn craton), and SEA09NC (Somerville et al., 2009—non-cratonic). Other models used in NSHA 2018 have been adapted from California (BEA14; Boore et al., 2014), Europe (CY08SWISS; Edwards et al. 2016 and CY14; Chiou and Youngs, 2014), and Central Eastern US (AB06; Atkinson and Boore, 2006). In previous versions of the NSHA, two other GMMs curated for Australia were used as candidates: Gea90SEA (Gaull et al., 1990—South-East Australia) and Gea90WA (Gaull et al., 1990—Western Australia. Gea90SEA and WA GMMs were not used in the 2012 (Leonard et al., 2014) or 2018 (Allen et al., 2018) NSHAs, as these GMMs were calibrated to local magnitude (M_L) and were developed based on the attenuation of macroseismic intensities. Future revisions of the NSHA are considering the adoption of the Next Generation Attenuation for Central and Eastern North America (NGA-East; Goulet et al., 2021) as a GMM candidate. For the purposes of this article, the “NGA-East” model refers to the weighted mean GMM from an ensemble of GMMs (Goulet et al., 2021) that we adjust from the reference site condition of V_S30 = 3000 m/s to a range of V_S30 values from 270 to 1100 m/s to account for the possible site conditions discussed herein. GMM constant terms and parameter uncertainties associated with V_S30 scaling represent sources of epistemic uncertainty (Ramos-Sepulveda et al., 2022) that we cannot resolve in this study.

Recently, Hoult et al. (2021a, 2021b) evaluated several GMMs in comparison to the ground motion recordings of the 2012 5.1 M_W Moe and 2021 Woods Point earthquakes in Victoria, Australia (Figure 1). For the Moe earthquake, there was a tendency for all GMMs to over-predict near-source (within a rupture distance of 40 km) ground motions; this was particularly the case for the mean NGA-East GMM (Goulet et al., 2017). Hoult et al. (2021b) concluded that AB06 and A12 perform well, at smaller epicentral distances (<70–100 km but >40 km) and particularly at shorter spectral periods. For the Woods Point earthquake, the absence of recorded ground motions at distances less than ∼60 km from the rupture precluded the evaluation of near-source GMM performance. Several GMMs (NGA-East, AB06, A12, Sea09) performed well at >60 km distances (Hoult et al., 2021b), with some variability in performance at different periods. At larger rupture distances (i.e. >100 km), GMMs that incorporate an attenuation “transition zone” due to the strong postcritical reflections from the Moho discontinuity (Allen et al., 2007; Atkinson, 2004) tended to provide better fits to the data at mid-to-longer spectral periods (e.g. AB06, A12, and NGA-East) (see Figure 2). This was also the case for the Moe earthquake (Hoult et al., 2021a). Hoult et al. (2021a, 2021b) emphasized the need to further evaluate GMMs to assist in logic tree weightings in PSHAs; we pursue this in this study.

Unreinforced masonry (URM) chimneys are among the most damage-prone components across any building class when subjected to earthquake ground motion (Krawinkler et al., 2012; Maison and McDonald, 2018; Moon et al., 2014). The natural period range of chimneys in 1–3 story buildings is 0.05–0.2 s. This period range falls within the typical high acceleration region of the design response spectrum (Klopp, 1996—see Chapter 3 therein); Maison and McDonald (2018) use a standard natural frequency of the host structure (i.e. house) of 6 Hz (0.16 s). Modeling of masonry chimney fragility curves can provide a prediction of seismic vulnerability and damage in the event of an earthquake as an expression of peak ground accelerations (PGAs) within a range of frequencies. This article presents a novel approach of using chimney damage states and their respective fragility curves to approximate earthquake shaking intensities using Bayesian techniques. Fragility curves define a cumulative distribution function with respect to a ground-motion intensity measure, commonly PGA. The intersection between the fragility curve and the output of a GMM defines the probability of damage (Heresi and Miranda, 2021; Maison and McDonald, 2018). In this study we use Bayesian modeling to evaluate the relative performance of GMMs based on their modeled PGA outputs relative to calculated PGA ranges from observations of chimney damage (or lack thereof) and associated chimney fragility curves.

Bayes’ theorem (Bayes, 1764; Joyce, 2003) states that one can calculate the posterior probability using a prior probability and likelihood function. If the likelihood of a calculated PGA (derived from a GMM) to cause damage to a chimney is known, the likelihood a GMM correctly predicts the observation can be calculated. Therefore, using Bayes theorem and residential chimney fragility functions, the likelihood that a GMM represents the expected damage outcome of a set of chimneys after an earthquake can be deduced.

This article aims to answer the following questions:

1. Using earthquake damage of chimneys and their respective fragility curves as a proxy for seismic ground motions at near-source locations, what is the relative performance of commonly used GMMs in terms of their ability to predict the chimney damage observations?

2. How do the relative weightings of GMMs established from a chimney analysis using a Bayesian approach compare with pre-Woods Point earthquake NSHA18 weightings of GMMs?

This article outlines a method to evaluate the relative performance of GMMs in an earthquake using proxy indictors (chimneys) to compensate for an absence of near-source instrumental seismological data. Two chimney fragility curve models (Model 1: Maison and McDonald, 2018; Model 2: Vaculik and Griffith, 2019) and median GMM PGA curves are used as inputs into a Bayesian analysis. The Bayesian approach is used to determine which GMM best predicts chimney damage/non-damage observations from the Woods Point earthquake.

At 9:15 am on 22 September 2021 (AEST), a moment magnitude (M_W) 5.9 earthquake occurred in southeast Australia within the Southeast Seismic Zone, approximately 130 km ENE of Melbourne’s CBD (Figure 1). This intraplate event was the largest earthquake in the state of Victoria since European record keeping began in the early 1800s (McCue, 2015). The epicentral region is sparsely populated; Woods Point (pop. 33), Jamieson (pop. 382), and Licola (pop. 11) (see Figure 4) are the three settlements nearest the epicenter (2021 Census). The mainshock occurred in the Victorian Highlands, approximately 13 km ENE from Woods Point. Four epicentral locations have been published (Table 1).

The location published by SRC is preferred as it uses additional data from industry and university seismic networks that are not used in the GA and USGS solutions. The north–south trending aftershock distribution (Figure 1) coupled with mainshock focal mechanism solutions suggests a steeply west dipping (83°; strike 172^o) (GA) to steeply east dipping (84°; strike 351^o) (USGS) left-lateral strike-slip source fault mechanism (Figure 1). SRC aftershock data were fit by linear interpolation in map view and cross-section to delineate an ∼8 km long, 9 km down-dip NNW-striking mainshock rupture plane with ∼85°E dip (Quigley et al., 2021; Quigley and La Greca, 2021). This rupture-source model is generally consistent with focal mechanism constraints and is used for the analyses herein. The nearest seismometer is a short period passive sensor, located 35 km from the epicenter at Thompson Reservoir (TOT) (Figure 1). Earthquake waveforms at TOT clipped under the Woods Point earthquake ground motions. Preliminary observations of ground motions (Hoult et al., 2021b) omit TOT data from analysis.

Geoscience Australia used four equally weighted GMMs (Atkinson and Boore, 2006 (AB06), Somerville et al., 2009 non-cratonic (Sea09NC), Allen, 2012 (A12), and Boore et al., 2014 (Bea14)) and a 0–30 m depth-time averaged shear-wave velocity model (V_S30) based on the Australian Seismic Site Conditions Map (McPherson, 2017) to produce PGA contour plots for the Woods Point earthquake (Figure 1) (e.g. Allen et al., 2019, 2023). These PGA contours are informed by the “ShakeMap” system in which GMMs, ground motion amplification (McPherson, 2007), and felt reports are combined to create a map of seismic intensity (Allen et al., 2019; Wald et al., 2010). Estimated PGAs within the epicentral region are ∼0.2 g. Calculated PGAs in Melbourne based on Geoscience Australia’s ShakeMap range from 0.02 to 0.05 g (Figure 1).)

Instrumentally recorded spectral accelerations (SAs) of the Woods Point earthquake are primarily within the range of the suite of median values predicted by GMMs at diverse periods and distances (Figure 2). The performance of GMMs at individual sites is challenging to evaluate because local V_S30 estimates are poorly constrained, and potential variations in V_S30 provide a source of epistemic uncertainty in evaluating GMM performance. Consequently, adherence to the assumed reference site condition of V_S30 = 760 m/s that is used to model the GMMs in Figure 2 and Hoult et al. (2021b) is uncertain; it is possible that inclusion of site-specific V_S30 data could affect the relative performance of GMMs (this is beyond the scope of this study but is a topic of ongoing research). SAs in excess of GMM estimates across the range of periods and source distances (Figure 2) could also reflect source effects associated with specific characteristics of the Woods Point earthquake (e.g. radiation effects, high source energy release, faulting mechanism, source geometry, path effects, and/or site effects such as amplification associated with low V_S30 host materials and basin effects) that may be distinct from the “average” observations from reference earthquakes used in GMM generation. Hoult et al. (2021b) show that short period (T < 0.1) SAs appear to be best matched by central and/or northeast American developed GMMs (NGA-East and AB06), while models developed for southeast Australia such as A12 “predict the attenuation, particularly at lower frequencies.” The absence of observational data at epicentral distances <60 km precludes comparison of GMM predictions against observations in the near-source region.

Modified Mercalli intensity (MMI) data derived from more than 43,000 “felt” reports range from MMI VIII to MMI I (Figure 3). Felt reports were submitted from epicentral distances ranging from <20 to >700 km (Figure 3a). Very few (<10) felt reports were submitted from the near-source region (<50 km) and there is large variability in MMI reports within “binned” epicentral distances (e.g. 40 km distance; Figure 3b), reflected by standard deviations (1σ) about the mean of >1 to 2 MMI increments. Mean MMIs generally decrease with increasing epicentral distance. Mean MMIs adhere most closely to the Australia-based (L15 AU) and eastern US-based (AWW14 CEUS) MMI versus distance intensity prediction models; AWW12 ATR (active crustal region) and AWW14 CA (California) models tend to underpredict observed MMIs (Figure 3b). At equivalent epicentral distances, MMIs tend to be higher northern and southern quadrants relative to eastern and western quadrants, although data are limited due to sparse population in the east. The dense urban population of Melbourne results in a sampling distribution biased into the western quadrant at 100–150 km epicentral distances. The dearth of near-source MMIs further illustrates the challenge of understanding near-source shaking intensities in this earthquake and highlights the need for consideration of other shaking intensity proxies.

A reconnaissance field survey of environmental and infrastructure damage in the epicentral area commenced ∼30 h after the earthquake (23–27 September 2021) (La Greca and Quigley, 2021; Quigley and La Greca, 2021). Subsequent field surveys were undertaken on 8–10 and 23–26 October. We identified 43 brick masonry and four stone chimneys (n = 47), including chimneys with observational evidence for damage (Figure 4). All chimneys were physically examined, precisely located, and photographed from multiple angles (Figure 5). Chimney width and heights were determined by brick counting on photographic images, using the average Australian brick size (76 mm high × 230 mm long × 110 mm wide), with exception of the four stone chimneys, where height was approximated using scaled photographs. Visual inspection of mortar in the field and in photographs enabled us to estimate mortar “quality” and “age.” This included up close investigation of the mortar, determining whether it was flaking and/or breaking apart. Mortar rankings from “Weak” to “Expected” were assigned to each chimney and entered as a data input in relevant fragility analyses (see fragility curve methods below). Mortar observations also enabled us to distinguish earthquake from pre-earthquake damage; for example, if a chimney had small cracking in the mortar and/or bricks but had moss growing within the crack, it was interpreted to be pre-seismic deformation. All chimney data are presented in Supplemental Appendices A and B. Of the total chimneys observed, five were determined to have collapsed in the earthquake or suffered extensive damage.

A schematic of the methodology is provided in Figure 6. The GMMs analyzed in this study are those used in the Australian NSHA18 and predecessor NSHAs, or are considered as candidate GMMs for future NSHAs (Table 2). The respective expert elicitation weights for NSHA18 GMMs are shown in Table 2.

The vulnerability of URM chimneys to earthquake ground motions was modeled using the fragility models of Maison and McDonald (2018) and Vaculik and Griffith (2019). The output of both models is a set of analytical fragility curves that express the probability of a chimney reaching a particular damage state as a function of ground motion intensity in terms of the PGA. Increasing chimney resilience shifts fragility curves toward higher PGA values. These curves in turn serve as the input into the Bayesian analysis in the subsequent portion of this article. The decision to consider two separate fragility models was made to improve the reliability of the Bayesian inference process given the inherent uncertainty in relating the expected damage states to ground motion intensity. PGA values are calculated using selected GMMs within the OpenQuake hazardlib software library at a distance equivalent to the source-to-site distance required by the GMM for each chimney on a range of site conditions (Pagani et al., 2014). These estimates were used to determine the probability of the chimney sustaining the degree of damage (or non-damage) that was observed. The Bayesian analysis represents the likelihood of each individual GMM to correctly predict the field-derived damage observations under the assumption that the fragility curve is representative of a chimney’s damage potential. In addition, the Bayesian analysis considers and incorporates two prior inputs. A uniform prior which assumes that all GMMs have an equal weight before entering the Bayesian analysis and an NSHA18 prior approach, where the NSHA18 logic tree weights are applied into the Bayesian analysis. PGA values output by GMMs in this process do not consider aleatory variability and assumes that the median GMM value is the PGA at the site of the chimney.

The fragility curves of Maison and McDonald (2018) were determined using a single-degree-of-freedom computer model (Maison and McDonald, 2018). Fragility curves incorporate the effects of various site parameters including chimney height above roof, masonry flexural tensile strength, chimney section dimensions, vertical steel reinforcement, and chimney-house anchorage strength. Damage functions for URM chimneys are expressed as a function of PGA, chimney height, and expected masonry tensile strength for the chimneys part of this study can be seen in Figure 7a and b.

A select set of chimneys in the Woods Point epicentral region did not meet the brick base templates used in Maison and McDonald (2018). Maison and McDonald (2018) suggest that the shortest measurement of either the length or width of the chimney dictates the damage function (Figure 8). Therefore, in this study for a chimney to be included it must meet only one measurement of the brick base outlined in Figure 7a and b.

Model 1 cannot be used for all chimneys due to non-adherence of chimney characteristics, including chimney brick type (stone chimneys were omitted), height, and brick base, to the parameters in the regressions (Figure 7); chimneys 6, 10, 11, 15, 20, 21, 22, 30, 31, and 41 were excluded from analysis. All chimneys are evaluated using Model 2 (see next section).

The damage function defines a PGA range of 50% probability of failure for varying tensile strengths and chimney brick base (Figure 7). To enable a range of values to be represented by a single value in fragility curve analyses, we extract the median PGA value of the 50% PGA range in the damage functions. We identify this as a model simplification as we do not consider the full range of the 50% chance of failure field. Following Maison and McDonald (2018), we use a beta value of 0.6 as per US FEMA P-58 Seismic Performance Assessment of Buildings, Methodology, and Implementation guidelines (Federal Emergency Management Agency (FEMA), 2018).

Model 2 follows the analytical two-step time-history analysis (THA) approach described in Vaculik and Griffith (2019). We first perform a THA on the parent building with excitation by the ground motion to compute the motion at the top of the building. This motion is then used to undertake a nonlinear THA of the chimney. The chimney’s force–displacement behavior was defined using a bilinear rule with a descending post-yield branch to represent rocking behavior (Vaculik and Griffith, 2017). Unlike the Maison and McDonald approach, this model ignores any bond strength and assumes that the chimney’s lateral load resistance is provided entirely from stabilization due to gravity. The force–displacement capacity of each chimney was constructed as a function of its geometry; that is, the height above the roof line and base width. In the case of rectangular chimneys (with unequal base widths), the shorter dimension was used. A factor of 0.9 was applied to the gross width of the chimney to account for deviation from idealized rigid behavior, for example, due to geometric imperfections and finite compressive strength.

Following the approach described in Vaculik and Griffith (2019), a set of five displacement-based damage levels were defined, ranging from D1 (first onset of cracking) to D5 (complete collapse). In order to align these with the observable damage levels in the field survey, these were condensed into three states: (1) no visible damage (damage < D2), (2) visibly damaged but not collapsed (damage ≥ D2 and ≤ D4), and (3) collapsed (damage > D4). Note that the onset of observable damage was set at D2 rather than D1, due to micro-cracking not being able to be visually assessed in the field.

This overall procedure was implemented within an incremental dynamic analysis (Vamvatsikos and Cornell, 2002) using a suite of code-compatible (Standards Australia, 2018) ground motions. The median PGAs to reach different damage states were computed, thus resulting in a standalone set of fragility curves for each chimney. Further detail regarding the overall approach can be found in Vaculik and Griffith (2019).

As with the first model, a beta value of 0.6 was adopted for the dispersion of the fragility curves consistent with FEMA guidelines.

Fragility curves were formulated in terms of the lognormal distribution, whose cumulative distribution function (CDF) can be expressed as:

where x is the PGA; Φ is the standard normal CDF operator; μ is the natural logarithm of the median PGA at each damage state; and β is the coefficient of variation (Lallemant et al., 2015), which was taken as 0.6 in both models.

Illustrative examples of fragility curves obtained using the respective models are shown in Figure 9 (see Supplemental Appendix A for all chimneys). These consider three chimneys:

• Chimney 17—a stocky (∼2 ft tall) chimney assumed to have weak bond (10 psi)

• Chimney 1—a medium-slenderness (∼5 ft tall) chimney assumed to have typical-strength bond (60 psi)

• Chimney 32—a slender (∼8 ft) tall chimney assumed to have normal-strength bond (60 psi)

It is seen that the solitary curves for Model 1 coincide roughly with damage state D5 curve in Model 2; and thus, the median PGA of the D4 curves in Model 2 (delineating the collapse state in the implementation throughout this article) are typically lower than the PGAs predicted by Model 1.

OpenQuake (https://platform.openquake.org/) hazardlib.gsim library was used to compute the expected PGAs at each chimney (Pagani et al., 2014). Chimney distances to the fault rupture were calculated and used as input into the GMMs called from the OpenQuake hazardlib.gsim library. Table 3 displays the fault geometry and inputs for GMM calculation.

The calculated PGA value from each GMM for each chimney can then be compared with the respective fragility curves. Figure 10 shows an example of five chimneys displaying that the intersection point of the expected PGA and the fragility curve determines probability values inputs for the Bayesian model.

The PGA–fragility curve intersection point represents the probability the chimney will exceed a specified damage state at that PGA for each given GMM for six selected chimneys. Depending on the damage state of the chimney resulting from the earthquake, the probability of the observed can be determined through either taking the intersection CDF (if chimney had damage), or through subtracting the CDF from 1 (if the chimney had no damage). M1 (Maison and McDonald, 2018) uses a binary damage state for collapse or no collapse (red line in Figure 10). M2 (Vaculik and Griffith, 2019) uses five damage states to model chimney damage (black lines in Figure 10). For chimneys that had damage but did not collapse, the CDFs for the damage range were subtracted from one another. These probabilities of observed values were used as inputs in the Bayesian approach.

Bayesian modeling is a statistical approach based on Bayes’ theorem using prior knowledge to update a statistical model (Van de Schoot et al., 2021). A prior distribution (prior knowledge) can be informed by new observational data to establish a new posterior probability. Bayes’ theorem states that one can calculate the posterior probability if the prior probability and the likelihood function is known. The approach used here calculates the relative probability of a GMM (A) correctly predicting “modeled” PGAs that are independently inferred from the chimney fragility curves and observed damage state of a chimney (B) as proxies for the “observed” PGAs at individual sites in the earthquake.

This can be expressed as:

where P(A|B) is the probability a GMM is correct given event B (the posterior belief), P(B|A) is the probability of the observed chimney damage state given the GMM output; the probability of observed. Table 4 outlines the probability of observed for each method and damage type.

P(B) is the summed probability of each chimney to be damaged as outlined by their fragility curve for all GMM and P(A) is the prior distribution of our GMM before the earthquake event. P(B|A) is determined through three different equations depending on the state of the chimney. As there are two methods of fragility curve generation, this analysis will examine the combination of these methods into a singular P(B|A) value.

There are two states of chimneys considered for Model 1; chimneys that have collapsed (P1), and chimneys that have not (P2). For chimneys that have collapsed, the fragility curve dictates the probability of the observed and therefore is represented as Equation 3:

For chimneys that have not collapsed, the probability is represented as the probability of not P1 (Equation 4):

The probability of observed for a GMM for fragility model 1 (P_observed(M1)) can be calculated as the product of P1 and P2:

Model 2 considers three damage states: undamaged, damaged, and collapse (Table 5). These are characterized by five damage levels (D-Levels) in which Table 5 outlines the relationship between the fragility curve D-Level and observed damage state of a respective chimney. Therefore, to calculate the P(B|A) for a GMM that incorporates three chimney damage states, three equations are derived. The undamaged state is represented in Equation 6. “No damage” is interpreted as not exceeding a D-Level of 2, and therefore the median value (⁠$μ$⁠) is defined by the D2 curves in Supplemental Appendix C for each chimney. Equation 7 is for chimneys that have sustained damage, but not total collapse. This is the probability the chimney would sustain a lower damage level, then subtracting the probability the chimney would fail from the higher damage level. This assumes a level of damage greater than a D-Level of 2 but not greater than 4. Equation 8 is for chimneys that have failed and therefore exceed a D-Level of 4. These three calculations state the probability of the observed event and can be seen in Supplemental Appendix D:

To account for the three different P(B|A) equations for the chimney damage states and therefore the probability of observed for M2 (P_observed(M2)), the P(B|A) for a given GMM is equal to the product of three equations. This can be expressed as follows:

This is completed for each GMM and provides a P(B|A) for both fragility curve methods.

P(B) is calculated from the sum of all P(B|A) values produced by all GMMs for each fragility curve model (Models 1 and 2). This results in an analysis that has two P(B|A) values determinant on the respective fragility curve model:

The analysis considers two iterations of priors. The first iteration assumes equal probability of GMM prior distributions. The GMM weights therefore reflect the likelihood of the respective GMMs based on only the Woods Point earthquake analysis (forthwith “uniform prior”). The second iteration takes the prior distribution from GMM logic tree weights determined from expert elicitation for NSHA18 (forthwith “NSHA18 prior”). The output of the NSHA18 prior model thus reflects an updated NSHA18 GMM weighting.

The above process was repeated using five different V_S30 parameters in the GMMs: 270, 400, 560, 760, and 1100 m/s. Multiple shear-wave velocities were considered and included in this study instead of a singular velocity due to not knowing the V_S30 at each chimney site. We therefore assume uniform V_S30 site conditions at the location of each chimney. To capture some of this uncertainty, the GMMs were adjusted using model-specific amplification factors based on varying V_S30 values. Some GMMs do not have amplification factors for varying V_S30. In these cases, ground motions calibrated to a reference V_S30 were adjusted using factors from other published studies. The amplification factors from Seyhan and Stewart (2014) were applied to A12 and SEA09NC. The NGA-East model was adjusted from the reference site condition of 3000 m/s using linear and nonlinear components of the amplification model developed by Stewart et al. (2020) and Hashash et al. (2020), respectively. The Stewart et al. (2020) model weights linear amplification contributions with large impedance contrasts relative to those with relatively gradual velocity gradients. Given typical site conditions for eastern Australian sites are more likely to possess more gradual velocity gradients (e.g. Collins et al., 2003), the full weight is placed on the gradual velocity gradient component of the Stewart et al. (2020) linear amplification model for adjusting the NGA-East model.

The Bayesian analysis calculates a P(A|B) for each of the two fragility curve models within the three Bayesian models (NSHA18 GMMs using uniform priors, NSHA18 GMMs using NSHA18 priors, and all GMMs using uniform priors). One source of epistemic uncertainty is associated with the performance of each fragility curve model and whether one model is more probable to accurately define the fragility of the sampled chimneys. This analysis assumes the fragility curve models are equally probable for four reasons: (1) chimney fragility curves are highly uncertain and there is insufficient information for a hierarchical analysis of both GMMs and fragility curves, (2) multiple chimneys could not be modeled by Model 1 and therefore result in a sample size discrepancy, (3) Model 2 is penalized for aiming to model chimney fragility damage at a higher resolution due to this fragility curve model having five damage states rather than only modeling collapse (Model 1), and (4) the Bayesian analysis is intended to evaluate GMMs and not chimney fragility curves. We hence average the P(A|B) value for both fragility curve methods for their respective GMM within each Bayesian model and therefore assume each chimney fragility model is equal.

Three Bayesian models were used to evaluate relative GMM model performance for ground motions produced by the Woods Point earthquake using damage to chimneys within the epicentral region. Results of the Bayesian analysis can be seen in Table 6 and can be visualized in Figure 11. The first two models (column 3 and 4 in Table 6; Figure 11a and b), isolate and evaluate only NSHA18 GMMs, the first using uniform priors, and the second using NSHA18 priors. The third model (column 5 in Table 6; Figure 11c) uses all GMMs mentioned in this study with uniform priors.

At the V_S30 values considered in this study, and considering only uniform priors, GMMs can be generally grouped by relative performance in the Woods Point earthquake. GMMs used in the NSHA18 for southeast Australia (A12, AB06, BEA14, CY08, CY14, and SEA09NC) perform relatively well compared with the other analyzed GMMs not used in NSHA18 (NGA-E, GEA90SEA, GEA90WA, SEA09YC) (Figure 11). This observation independently supports the selection of GMMs for NSHA18 (Griffin et al., 2018).

Based on the predominance of sedimentary bedrock in the study area, and on the topographic slope as reflected in the seismic site condition map for Australia (McPherson, 2017; Wald and Allen, 2007), the most likely representative site class for the study region is B to B/C (V_S30 of 760–1100 m/s). In this V_S30 range, A12, AB06, and CY08SWISS are the best performing GMMs, both with uniform and NSHA18 Bayesian priors (Figure 11, Table 6). However, several of the chimney sites are on alluvial flood plains of Quaternary gravel (e.g. Licola, Jamieson) or on sloped sites (e.g. Woods Point). These spatial variations are predominantly smaller than the spatial resolution used to determine V_S30 pixel values for the site class condition map. The thickness of low-velocity unconsolidated sediments at each site is unknown (probably <5–10 m, but not measured) and the potential effect of slope on site amplification is an additional epistemic uncertainty. To partially address this uncertainty, we also considered relative performance of GMMs at lower V_S30 values (270–560 m/s). CY08SWISS, CY14, and AB06 are the best performing GMMs in this V_S30 range with uniform and NSHA18 priors (Figure 11, Table 6). Further study could include estimation of V_S30 at each individual site and consideration of nonuniform V_S30 parameters within a single model; we acknowledge the simplification of uniform V_S30 parameters within individual models as a limitation of our study.

A comparison of NSHA18 GMM performance in the Woods Point earthquake (uniform priors) with the pre-earthquake NSHA18 GMM weightings (NSHA18 priors) is useful to consider how the outputs of the Bayesian approach (which could inform GMM weightings in future NSHAs) relate to the NSHA18 weightings (Table 7). The percent residual is a metric for this difference; a positive residual reflects an increase in the relative performance (weighting) of a GMM through incorporation of the Woods Point performance with the NSHA18 priors, and a negative residual indicates a reduction in relative performance. In the most likely V_S30 range, some GMM percent residuals are dependent on which V_S30 value is used (e.g. A12, CY14), while for others, weightings either increase (e.g. AB06, CY08SWISS) or decrease (BEA14, SEA09NC).

There are substantive uncertainties involving the implications of the refined GMM weightings using the Woods Point earthquake (Table 7) to weighting GMMs in future NSHAs. These include (1) the GMM relative performance for the Woods Point earthquake is determined only for this earthquake, and thus the applicability to GMM relative performances in future earthquakes is uncertain; (2) the study is confined to the near-source region, and it is therefore unknown whether the weightings are relevant at distances beyond the near-source; (3) the study sites are not uniformly spatially distributed with respect to the earthquake source and thus could be susceptible to spatial anisotropies in ground motions and associated sampling biases; (4) the faulting kinematics of the Woods Point earthquake are distinct from prior events that informed NSHA18 GMM weightings, and thus diverse faulting mechanisms in future earthquakes could result in different weightings than those presented here; (5) possible spatial variations in V_S30 among sites are not considered, and thus GMM relative performance rankings could vary if site-to-site V_S30 variations were applied; and (6) GMM performance is evaluated using non-instrumental proxies and thus translation to specific instrumental measures including PGA and SA is challenging. This issue of using a single event to evaluate GMMs is further complicated by the inter-event variability of Woods Point earthquake. Observed ground motions at different periods may be substantially higher or lower than median predicted ground motions (Baker and Cornell, 2005). GMM—instrument analysis in Figure 2 and Hoult et al. (2021b) suggest that at distances of 50–150 km, the Woods Point earthquake is a positive epsilon event (i.e. has specific period ranges where SAs exceed median spectral response curves; Baker and Cornell, 2005). If this event is a positive epsilon event, our study would preference GMMs that estimate higher PGAs and ultimately affect the final weightings. A comprehensive approach would be to repeat this analysis over multiple earthquakes, as randomly chosen events should ultimately aggregate to an epsilon averaging zero (Baker and Cornell, 2005). A final uncertainty is the fragility curves themselves. Different methodologies, chimney parameter inputs (e.g. masonry tensile strength), and beta values will impact the results and interpretations of this study.

Acknowledging these uncertainties, consideration could be given to down-weighting GMMs that under-performed in the Woods Point earthquake and increasing the relative weights of GMMs with increased performance for future NSHAs; however, we do not suggest that absolute weights assigned in any future GMM logic trees should adhere resolutely to the Bayesian output weightings in Table 7. Given the epistemic uncertainties, aleatoric variability, and potential for sampling biases, we cannot advocate for exclusion of any of the NSHA18 GMMs or NGA-East from consideration in future NSHAs on the basis of their relative performance in the Woods Point earthquake alone.

In terms of specific GMM performance, A12, AB06, and CY08SWISS are the GMMs that best match the chimney damage observations. At low V_S30 values (270 and 400 m/s), the clear statistical preference is CY08SWISS, presumably because V_S30 scaling in other GMMs overestimated PGA relative to the PGA inferred from observations of chimney damage. This is attributed to CY08SWISS predicting relatively low PGAs at low V_S30s. At high V_S30s, A12 is the best performing GMM. The A12 model is curated for southeastern Australia and uses a stochastic finite-fault simulation technique (e.g. Motazedian and Atkinson, 2005) involving the use of reinterpreted source and attenuation parameters for small to moderate magnitude southeast Australian earthquakes. Similarly, the AB06 GMM uses a stochastic finite-fault simulation technique using earthquakes from the ENA region. A12 and AB06 models are both considered to perform well due to similar simulation techniques and the use of similar geometrical attenuation coefficients in the near-source region (Allen, 2012; Allen and Atkinson, 2007; Atkinson, 2004; Hoult et al., 2021a). Allen and Atkinson (2007) concluded that there is no significant difference in source characteristics of ENA and SEA earthquakes. This was used for justification of the inclusion of AB06 in the NSHA18 GMM suite. CY08SWISS also performed well. The CY08SWISS model was adjusted for use within the 2015 Swiss Seismic Hazard map (Edwards et al., 2016) and this variation of the model has been used here. The geology of deformed and thrusted Mesozoic sediment overlying crystalline basement in Switzerland (Pfiffner, 2021) is crudely geologically analogous with the deformed Silurian to early Middle Devonian Victorian Highlands of Australia (Fergusson et al., 1986) and thus, similarities in seismic attenuation characteristics between the two regions could be expected. The SEA09NC model performs relatively poorly compared with other NSHA18 models across the range of V_S30, but outperforms non-NSHA18 models at all V_S30 values excluding 1100 m/s (Figure 11). If the performance in the Woods Point earthquake is borne out in other earthquakes, future consideration of NSHA logic tree weightings could consider downweighing or exclusion of the SEA09NC model for this part of Australia. The Yilgarn craton (SEA09YC) version of this GMM performs poorly and on average predicts the second highest expected PGA behind NGA-East, which should have resulted in more chimney failure. The CY14 and BEA14 models represent the California region of the Western US. These models perform relatively well in the 560–760 m/s velocity range. However, the performance of these models—particularly BEA14—decreases at the 1100 m/s range. The CY14 model uses Z_1.0 and Z_2.5 inputs intended to model basin effects. However, these variables were not measured at each chimney site. Rather, standard generic proxies are used. This omits Z_1.0 and Z_2.5 as contributing factors into the increased likelihood of this model within low velocity ranges (400 m/s). This highlights CY14 as a well-performing GMM without consideration of basin effects. The Gaull et al. (1990) models (SEA and WA) both performed poorly. These models were not included in NSHA18 due to poor model performance against ground motion records and challenges with the conversion of local magnitude to moment magnitude. Their poor performance in the Woods Point earthquake further justifies their exclusion from future NSHAs. The NGA-East model is the among the worst performing models for the near-source region. We find that NGA-East likely overestimates ground motions within the epicentral region for the Woods Point earthquake.

Considering all of these aspects, one hypothesis that we are continuing to test is that ground accelerations in the near-field region of the Woods Point earthquake in the frequency range relevant to chimney damage were lower than anticipated relative to “median” ground accelerations for an earthquake of this magnitude. Other environmental features may shed light on this hypothesis (La Greca and Quigley, 2021; Quigley and La Greca, 2021).

We developed a Bayesian approach to evaluate the relative performance of GMMs in the near-source (0–50 km) region of the moderate magnitude (5.9 M_W) intraplate 2021 Woods Point earthquake using independent observations of chimney damage and chimney fragility curves. For each chimney, the probability of a given GMM producing PGAs consistent with chimney damage states was calculated across a range of V_S30 values. At the most likely V_S30 (760 m/s), the best performing models are AB06, A12, and CY08SWISS. GMMs that were preferentially selected for utility in the Australian National Seismic Hazard Model (NSHA18) prior to the Woods Point earthquake outperform other GMMs, providing additional, independent support for their selection. The recently developed NGA-East GMM performs relatively well in the more distal region (e.g. >50 km) but is among the poorest performing GMMs in the near-source region across the range of V_S30. This could reflect unsuitability of the V_S30—adjusted NGA-East GMM relative to other GMMs to predict near-field ground motions for the Woods Point earthquake and possibly other earthquakes in southeast Australia. Further study of the parameters and uncertainties in NGA-East GMMs and site amplification models may provide further insights into this observation. Alternatively (or in addition), the relative performance of GMMs could reflect epistemic uncertainties in our proxy approach (e.g. methodological deficiencies and/or parametric uncertainties in using chimneys to estimate ground motions). A final consideration is that ground accelerations in the near-field region of the Woods Point earthquake in the frequency range relevant to chimney damage were lower than anticipated relative to “median” ground accelerations for an earthquake of this magnitude. The Woods Point earthquake is the largest instrumentally recorded earthquake in this region and the first moderate magnitude (M_W > 5) strike-slip earthquake in Australia; as such it provides valuable information to enrich the seismic catalog and consider the relative performance of GMMs in predicting ground shaking intensity distributions here and analogous regions globally.

We want to thank the Seismology Research Centre for sharing their seismic data, which greatly contributed to this study. We’re also grateful to AuScope (www.auscope.org.au) for providing funding for earthquake reconnaissance that allowed us to collect the data used in this manuscript.