A Bayesian network‐based probabilistic framework for updating aftershock risk of bridges

The evaluation of a bridge's structural damage state following a seismic event and the decision on whether or not to open it to traffic under the threat of aftershocks (ASs) can significantly benefit from information about the mainshock (MS) earthquake's intensity at the site, the bridge's structural response, and the resulting damage experienced by critical structural components. This paper illustrates a Bayesian network (BN)‐based probabilistic framework for updating the AS risk of bridges, allowing integration of such information to reduce the uncertainty in evaluating the risk of bridge failure. Specifically, a BN is developed for describing the probabilistic relationship among various random variables (e.g., earthquake‐induced ground‐motion intensity, bridge response parameters, seismic damage, etc.) involved in the seismic damage assessment. This configuration allows users to leverage data observations from seismic stations, structural health monitoring (SHM) sensors and visual inspections (VIs). The framework is applied to a hypothetical bridge in Central Italy exposed to earthquake sequences. The uncertainty reduction in the estimate of the AS damage risk is evaluated by utilising various sources of information. It is shown that the information from accelerometers and VIs can significantly impact bridge damage estimates, thus affecting decision‐making under the threat of future ASs.

The accurate evaluation of the state of infrastructural assets in the aftermath of major earthquakes is of paramount importance for optimal emergency management decisions and ensuring their safe use under the threat of ASs. It is also essential to correctly allocate and prioritise resources to minimise further casualties and business interruptions and speed up recovery from disruption (e.g., 2 ).
Many approaches and frameworks have been proposed in recent decades to address the post-earthquake functionality of bridges and their AS risk, aiming to support the decision on opening the bridge to traffic. Mackie and Stojadinović, 3 among others, developed an approach based on the performance-based earthquake engineering (PBEE) framework by the Pacific Earthquake Engineering Research (PEER) Center to evaluate the loss in load-carrying and traffic-carrying capacity of bridges following an MS. However, their approach disregarded the increased risk of damage due to ASs. On the other hand, Franchin and Pinto 4 proposed a method for deciding whether or not to allow traffic on a bridge based on the survival probability of an MS-damaged bridge under AS hazard. However, the problem of the uncertain knowledge of the bridge state following an MS was not addressed in their study. Instead, they assumed that the bridge damage is both analytically and visually detectable and that a near real-time estimation of the MS source and site parameters is available (e.g., from regional seismic stations/networks). The accurate evaluation of the level of damage suffered by a structure following an earthquake is critical for a correct estimation of its AS risk (e.g., 5,6 ). Obviously, in the absence of direct evaluations of damage (e.g., via field inspections), knowledge of the earthquake-induced ground-motion intensity experienced by the structure is essential for inferring the level of damage sustained by it. The simplest approach for evaluating ground-motion intensity at a given target site involves using a ground motion model (GMM), estimating the amplitude of a ground-motion intensity measure (IM) given the information on the location and magnitude of the earthquake (and eventually other site-specific features). However, the dispersion associated with GMMs is generally large, thus leading to highly uncertain estimates of the ground-motion IM. This uncertainty could be reduced by using information from seismic stations near the site (e.g., [7][8][9][10][11] ). However, the correlation between intensities at two points reduces significantly with the distance between them; hence, the IM estimate is generally still quite uncertain unless a station is located very close to the site. In any case, a fragility model (i.e., probability of damage as a function of a hazard intensity measure) is needed to estimate the damage level given the uncertainty of the IM (and modelling uncertainties). This results in uncertainty propagation in damage and loss estimates, leading to further uncertainty in assessing structural/nonstructural risk under future earthquakes.
One way to improve the knowledge of the actual state of a structure of interest is by exploiting observations from structural health monitoring (SHM) systems (e.g., 12,13 ). Most existing SHM methodologies rely on vibration measurements through accelerometers to detect changes in the dynamic properties of the system that can be attributed to structural damage (e.g., 14 ). Alternative sensors have also been proposed for the dynamic identification and damage detection (e.g., Global Positioning System (GPS) receivers, cameras, etc.) (e.g., [15][16][17][18]. Visual inspections (VIs) after earthquakes help determine the damage's severity, though they only partly capture effects of the seismic action (e.g., surface cracks in concrete, residual displacements [RDs]). Moreover, conclusions on the level of damage associated with these effects can be subjective, prone to human error and unreliable.
Further studies have proposed approaches for improving bridge damage estimates using information from sensors and VIs. For example, Alessandri et al. 19 proposed a method for evaluating post-earthquake bridge operability based on a rational combination of information derived from numerical analyses and in situ inspections. A Bayesian approach was developed to update the MS damage probability and thus the AS risk using the observation made by VIs. The study showed that in situ inspections could drastically affect and modify the damage estimates for a structure, helping to decide whether to open traffic on bridges.
Limongelli et al. 13 developed a framework for quantifying the value of SHM sensor information for post-earthquake emergency management of bridges, discussing alternative monitoring strategies. A similar framework was applied by Giordano and Limongelli 20 to evaluate the value of information of SHM sensors for bridge risk management, considering AS hazard. The first natural frequency of the structure was considered as damage sensitive feature in the study. Tubaldi et al. 21 developed a Bayesian network (BN)-based approach for evaluating and comparing the contribution of alternative data sources such as free-field seismic stations, GPS receivers, and structure-mounted accelerometers for bridge seismic risk assessment purposes. The comparison was conducted using two alternative measures that quantify the added value of information from the observations, based on pre-posterior variance and relative entropy reduction concepts. It was shown that the information from an accelerometer mounted on a bridge deck is superior in terms of uncertainty reduction to that provided by seismic stations located not very close to the site or by GPS receivers.
The assessment of the AS risk for bridges -and thus decision making on post-earthquake emergency management operations -may significantly benefit from the synergy of the various approaches outlined above and from the fusion of heterogeneous pieces of available information. The BN framework developed by Tubaldi et al. 21 only considers an MS scenario; in this study, the framework is further extended to include ASs in the risk assessment of bridge structures. Furthermore, for the first time, risk updating is carried out by combining valuable information from seismic stations, SHM sensors and VIs. The framework extension requires defining an AS hazard model that describes the frequency and intensity of ASs following an MS and developing a model that quantifies the damage accumulation under multiple ASs, given the damage experienced in the MS and the IMs of the ASs. The framework is applied to evaluate the AS risk of a case-study two-span bridge model, which is assumed to be located in Central Italy. This area was exposed to a series of ASs during the Central Italy earthquake of 2016-2017 (e.g., 22,23 ). Section 2 illustrates the various models involved in the proposed BN. Section 3 presents the BN developed to describe the relationship between the various parameters involved in seismic damage assessment and updating of these parameters based on additional available information from different sources. Finally, Section 4 illustrates the implementation of the method on the two-span bridge as an illustrative case study. This is followed by a discussion of the results and some concluding remarks in Section 5.

MODELS FOR AFTERSHOCK RISK ASSESSMENT
This section describes the various models required to develop the BN for AS risk assessment, along with the parameters involved and the observations required for updating them.

Mainshock analysis
A GMM is required to estimate the probability distribution of a ground-motion intensity measure for the MS experienced by the considered bridge at the ith site of interest, , given the following variables that are assumed to be known at the end of the earthquake: moment magnitude of the MS earthquake (M MS ); a measure of the source-to-site distance R MSi , for instance, assuming a point-source event, for simplicity; other parameters, collected in the vector s i , characterising the fault (such as those describing the faulting mechanism, fault geometry, depth to the top of the rupture) and site properties. In general, a GMM is characterised by the following form (e.g., 24 ): where ( , , ) is a function describing the lognormal mean of given M MS , R MSi and s i ; η the inter-event (or between-event) error term and is the intra-event (or within-event) error term.
The inter-and intra-event variabilities of a GMM describe the multi-level variabilities inherited in the ground motion recordings that arise from multiple station recordings for multiple earthquake events. Specifically, the inter-event error term describes the systematic variability in the ground motions throughout the region produced by different earthquakes of the same magnitude and rupture features. The intra-event error describes the variability in ground-motion intensity at various sites of the same soil classification and distance from the source during a single earthquake (e.g., 25 ). Thus, following Park et al. 26 and Crowley et al., 27 the same inter-event variability is applied to all sites of interest for a given earthquake scenario. In contrast, the intra-event variability is represented by a spatially correlated Gaussian random field. This can be built based on the intra-event error terms and the correlation coefficient ρ ij between the ground-motion parameters at two sites i and j for i,j = 1,2,..N sites , where N sites is the number of sites of interest. The corresponding covariance matrix of the ground motion IM field has the following form: where σ η and σ ξ represent the standard deviations of the inter-and intra-event error terms provided by the GMM, respectively. Further details about this representation of the ground-motion field can be found in Gehl et al. 8 The field observations of the ground-motion parameters at seismic stations can be used as evidence to update the prior estimates of the IM at the site of interest. The spatial correlation structure between the IMs at the monitored points and the target site plays a significant role in propagating the information from observed IMs to unobserved ones (e.g., 28 ).

Aftershock hazard analysis
The magnitude and location of an MS event can be used to obtain realistic AS scenarios. Several alternate methods available in the literature can be used to compute the seismicity/source parameters of ASs from the parameters of MSs. For example, Reasenberg and Jones proposed a stochastic Bayesian parametric model that allows estimating the probabilities of occurrence for ASs and larger MSs during intervals following an MS. The estimate of the model parameters is obtained with Bayesian statistics using the ongoing AS sequence and a suite of historic California AS sequences. Currently, one of the most popular models to describe the time-dependent seismicity of a region is the space-time Epidemic-Type Aftershock Sequence (ETAS) model (e.g., 30,31 ). The model's premise is that all earthquake events tend to trigger ASs with a stationary magnitude-dependent relation between them. The model considers earthquake seismicity from ASs (caused by internal stress adjustments in the seismogenic system) and background earthquakes (caused by forces within the plate tectonics, fluid/magma intrusion, slow slip, etc.). The ETAS model is a point process representing the occurrence of earthquakes for a given magnitude threshold over a given temporal and spatial range. Due to the computational complexity of the ETAS model, particularly in terms of the location of ASs, simpler approaches such as the branching aftershock sequence (BASS) model (e.g., 32 ) are available that probabilistically compute consistent source-to-site distances between the MSs and ASs. Apart from the models that relate the earthquake sequences using the causal event parameters, some methods propose the casual relations based on the intensity measures of the ground motions. For example, Fayaz et al. 33 proposed a novel approach that uses time-series modelling concepts to temporally correlate the Arias intensity (used as a proxy for the energy content of ground motions) of various earthquakes sequences. Generally, for the purposes of structural analysis and risk analysis, the above-mentioned procedures are also used to select consistent MS-AS ground-motion scenarios and records. Recently, Goda and Taylor 34 proposed an AS record selection procedure that can be used to simulate time-series data for MS-AS sequences. The study used the generalised Omori's law (e.g., 35 ) as the basis of generating artificial sequences and simulated a series of events with their time and magnitude stamps. For simplicity, they used the same source-to-site distance R for MS and AS events. Alessandri et al. 19 combined AS hazard analysis using Omori's law and Latin Hypercube Sampling technique for selection of MS-AS ground motions to conduct an AS risk assessment of bridge structures. Furthermore, a simplistic approach for simulating MS-AS sequences is to use Båth's law (e.g., 36 ) that postulates that the largest AS for a given MS is on average 1.2 moment-magnitude units lower than the MS.
In this study, for simplicity, the general procedure adopted by Papadopoulos et al. 37 is utilised to compute the number of ASs, their moment magnitudes ( ), time between the MS and ASs (Δ ), and the distance between the epicentre of the MS and that of the AS ( epi ). The sequences are generated by means of the triggering component of the ETAS model. The model is based on three uncoupled probability distributions that model (a) the direct offspring productivity; (b) the spatial; and (c) the temporal distribution of the triggered earthquakes, as well as a magnitude distribution of earthquakes derived from the Gutenberg-Richter (GR) law. Specifically, the simulation of MS-consistent AS scenarios in this study uses Equations (3)- (6). First, the number of direct offspring events from an MS event with a magnitude is sampled from a Poisson distribution with mean ( ) given in Equation (3). Then, for each of the offspring events, the distance epi between the epicentres of the two events and the inter-arrival time Δ between the parent and offspring event are simulated using Equations (4) and (5), respectively, where u t and u r are uniformly distributed random variables over the range [0,1], and A, a, p, c, D, q, γ are constant parameters representing the general spatial and temporal distributions of AS events and are estimated mainly through maximum likelihood estimation (e.g., 31,38 ). Similarly, for each offspring AS, the corresponding magnitude is sampled from the GR distribution using Equation (6), where u m is a uniformly distributed random variable over the range [0,1], b a constant parameter that represents the regional level of seismicity and min is the minimum considered magnitude. Equation (6) refers to the case of untruncated GR distributions with only a lower bound. After the direct offspring ASs (i.e., events triggered by the MS) are defined, the second generation of offspring events (triggered by direct ASs) can be sampled by repeating the procedure using the first-generation offspring as parent events. The second-generation offspring events can be again used as seeds for a third set and so forth until the sequence eventually dies out (zero offspring are sampled) or there are no more seed earthquakes within the time of interest. A reasonable time span can be used to dictate the time of interest within which the majority of ASs is expected to occur, and bridges are not likely to be repaired (e.g., 1-3 years). In this study, whose focus is on the rapid response to ln(EDPi) Illustration of the bilinear regression model used in this study earthquakes, a shorter time span of a few weeks might be considered. Earthquakes simulated outside of this time interval are usually discarded (e.g., 37 ).
The same GMM employed for the MS can be used to relate the IM of the ASs, collected in the vector IM AS , to the AS source to site distance (R AS ) as well as the magnitude of the ASs (M AS ). However, the framework can be easily updated to include AS specific GMMs, based on studies like Lee et al.; 39 similarly, the correlations between and within the IM estimates of MS and AS ground motions can be easily included in the proposed approach. Yet, the use of average spectral acceleration over a period range as the considered IM reduces significant differences. Furthermore, there is a lack of usable AS recordings in the available databases to develop AS-specific GMMs and correlation models (e.g., 37 ). Hence, most similar studies (e.g., 40 ) have reverted to the use of the same GMMs for both MS and AS.

Mainshock response and damage assessment
Given the MS ground-shaking intensity, a statistical model is required to describe the joint probability distribution of the engineering demand parameters (EDPs) of interest for increasing intensity levels (e.g., 41 ). Furthermore, modelling the correlation structure between the various parameters of interest is very important as it is the basis for updating the probabilistic distribution of one EDP (e.g., floor acceleration in a building) given the observation of another (e.g., storey drift). Alternative approaches can be employed to develop a joint probabilistic seismic demand model (PSDM), such as multistripe analysis (e.g., 42 ), incremental dynamic analysis (e.g., 43 ) or cloud analysis (e.g., 44 ). Cloud analysis is adopted in this study, given its computational efficiency and resulting accuracy. For this purpose, the computational model of the considered structure is analysed under a set of ground-motion records of different IM levels. The response parameters (EDP i , for i = 1, 2, .., N EDP ) are used as target variables for a regression model. In particular, a bilinear model is considered in this study (e.g., 45,46 ) since it allows a better description of the trend of the structural response with the ground-motion intensity. The model for the generic ith EDP is given in Equation (7) (see Figure 1):  damage index classification

Level
Calculated damage index Degree of damage where a 1 is the intercept of the first segment, b i for i = 1, 2 are the slopes of the two segments, IM* is the breakpoint IM, which is defined as the point of intersection of the two segments. The step function H(⋅) controls which of the two segments must be considered (i.e., H = 0 for IM ≤ IM*, and H = 1 for IM > IM*). The probability distribution of each EDP is also described by the values of two random variables (i.e., the error terms 1 (∼ (0, 1 )) and 2 (∼ (0, 2 ))), which are characterised by a normal distribution with zero mean and standard deviations 1 and 2 . Moreover, to define a joint probability density function (PDF) for the various EDPs, a covariance matrix must be assigned, which has the same form as that of Equation (2). For this purpose, different correlation coefficients must be estimated for the two conditions corresponding to IM ≤ IM* and IM > IM*, thus, leading to two correlation matrices. The EDPs considered in this study are the peak transient displacement (TD), the RD and the peak absolute acceleration (PA) at the bridge deck level. The first two parameters may be measured using GPS receivers, laser vibrometer, transducers and cameras, whereas PA is measured using accelerometers.
In addition to the above response parameters, two more EDPs are needed to define the VI outcomes, and , respectively, denoting the maximum strain of concrete cover in compression and tension experienced during the MS. For simplicity, it is assumed that a VI provides information only on whether concrete cracking and crushing have occurred or not at the base of bridge piers. These two events are defined by and exceeding their respective limits̄c c and̄c t , which are in general also random variables. The state of the bridge is described in this study by a global damage index (D), which is an extension of the Park and Ang 47 damage index to the case of biaxial loading. The Park and Ang 47 damage index is expressed by: where is the maximum displacement of the structural member, the ultimate curvature displacement, ℎ the dissipated hysteretic energy and is the yielding strength of the structural member; is a dimensionless empirical factor describing how hysteretic energy contributes to damage compared to the displacement demand. Experimental values of are in the range between −0.3 and 1.2 (e.g., 48 ). Table 1 summarises the classification of damage levels based on Park et al., 49 for establishing a correspondence between the calculated values of D and the observed empirical damage.
Since the structure examined is subjected to biaxial loading, a biaxial damage index is employed to define the bridge state. Rodrigues et al. 50 carried out an experimental campaign on 24 reinforced concrete (RC) columns tested under bidimensional earthquake conditions and presented seven expressions for the evaluation of D. The present study considers the following one: In Equations (10) and (11), and refer to the damage indices calculated for each independent direction, and are the cumulative dissipated energy for each independent direction, max, and max, denote the maximum displacements, ult, and ult, represent the ultimate displacements and yield, and yield, denote the yielding strength in each direction. To exploit the information from sensors and the outcomes of VIs for updating the knowledge of the bridge damage, the joint probabilistic model must include both the EDPs of interest and the damage index.

Aftershock damage assessment
The process describing the evolution of damage in a bridge under repeated events such as an MS-AS sequence is a statedependent process, that is, the increment of damage that characterises a given shock depends on the history of damage accumulated during the previous shocks. Previous studies (e.g., 51 ) considered the Markovian assumption, that is, the damage evolution under a given ground motion, conditional to the features of the ground motion, depends only on the state of the structure at the time of the shock and not on all of its damage history. Ghosh et al. 52 demonstrated in their study that earthquake damage accumulation can be treated as a homogeneous Markov process. Introducing this assumption, the damage index at the end of the nth event or ground motion can be expressed as a function of the intensity of the event and of the damage at the end of the (n-1)th event or ground motion. Specifically, the multi-variate model developed by Ghosh et al. 52 is employed in this study, which has the following form: where is the ground motion intensity of the nth event, −1 is the damage index after the (n-1)th event; , , , and are regression coefficients, and (∼ (0, )) is a random variable normally distributed with zero mean and lognormal standard deviation . This model can be fitted to empirical data from analyses of the structure under a series of MS-AS sequences, where = for n = 1. It is noteworthy that Equation (12) may return values of lower than −1 due to the nature of the regression model and . One way to overcome this physical inconsistency (i.e., damage can only increase) is to postulate that > −1 : This choice is expected to provide more physics-based results than the model originally developed by Ghosh et al. 52 It is worth mentioning that under the assumption of a homogeneous Markov process, the coefficients of the model in Equations (12) and (13) are independent of n, that is, the probability of moving from a state of damage to another under a given ground motion does not depend on the number of ground motion in the sequence (e.g., the same model can be used for the second earthquake and the third earthquake in a sequence).

BAYESIAN FRAMEWORK
This section presents the Bayesian framework developed for the assessment of the AS risk of a bridge and for its updating using observations from seismic stations, structural monitoring sensors and VIs.

Bayesian network
This subsection illustrates the BN developed to describe the probabilistic relationship between the parameters specified in the previous section, perform predictive analysis and update these parameters based on additional information from different observations (see Figure 2). The magnitude of the point-source MS earthquake (M MS ) and location between the source and the site as well as the source and the seismic stations (collected in R MS ) are assumed to be known. Various types of information are considered to be available for updating the probabilistic relationships of the variables in the network: on-site seismometers located close to the site of the structure, providing information on IM levels; GPS data, updating the knowledge of the RD; accelerometer data, updating the knowledge of the PA in the bridge deck; and the outcome of a VI, updating the knowledge of and . In the last case, the values of and are updated only if concrete cracking F I G U R E 2 Bayesian network illustrating the relationships between the parameters involved in the damage assessment (observed quantities are indicated by thick lines, parent nodes filled with grey) and/or crushing are observed, as this means that the deformations in the most critical fibre of the base section have likely exceeded the limit threshold corresponding to tensile cracking and/or crushing. Since the BN is static, it does not account for the temporal evolution of the system between the end of the MS and the occurrence of the ASs. Thus, information from sensors and VIs gained at different times after the occurrence of an earthquake can potentially be merged together. The nodes of the BN represent random variables, each characterised by a PDF. Nodes are related to their parent and child variables through edges that state conditional dependencies between variables (i.e., use of conditional probability distributions). Nodes with no parents are termed root nodes, associated with marginal probability distributions. Two forms of probabilistic inference can be carried out in BNs: predictive analysis that is based on evidence (i.e., information that the node is in a particular state) on root nodes, and diagnostic analysis, also called Bayesian learning, where observations enter into the BN through the child nodes. When evidence enters the BN, it is spread inside the network, thereby updating the probability distributions of the variables through one of the two forms of inference mentioned above.
The earthquake-induced ground shaking for the MS event is modelled by the deterministic root nodes M MS and R MS . For demonstration purposes, two seismic stations (represented by IM MS2 and IM MS3 ) are assumed to be in the vicinity of the bridge site (represented by IM MS1 ).
Following Gehl et al., 8 the inter-event variability is modelled by the root node W, which is parent to the three IMs of interest and follows a standard normal distribution. The intra-event variability is modelled via three root nodes, U j , for j = 1,2,3, also following a standard normal distribution. The following relation expresses the joint conditional distribution of the IMs given W and U i : where is the median value of the IM at the ith site, and ln ( , ) is the lognormal mean, which is a function of and (see also Equation 1), and are the lognormal standard deviations describing the intra-and inter-event variability, respectively, is a term of the lower triangular matrix obtained through a Cholesky factorisation of , which is the spatial correlation matrix expressing the correlation between the IMs at the various sites.
A similar approach is used for the PSDM describing the conditional distribution of the EDPs and damage index D MS given the IM at the site, IM MS . However, in this case, a bilinear regression model is employed, and thus two different error variables and correlation matrixes have to be considered, one for < * and the other for > * . In addition, three additional root nodes, denoted as e VD , e GPS and e PA , are introduced. They are used to describe the errors in the visual damage estimation and the measurement errors of the observations obtained with GPS and accelerometers. These error variables are assumed to be zero-mean normally distributed variables. They implicitly represent the accuracy and reliability of the observations from sensors and the VI: for example, if the standard deviation of e VD is large, the information from VIs is not very reliable. Thus, its contribution in the BN updating will be almost negligible.
The IMs of the forecasted AS events, collected in IM AS , are computed by using the same GMM used for the MS, based on the knowledge of the distances between the AS source and the site, collected in the vector R AS , and the magnitude of the AS sequence M AS . Finally, the model described in Equation (12) is employed to describe the relationship between the damage incurred from the MS D MS and the damage incurred from the ASs D AS , given the IMs of the forecasted AS events collected in IM AS .

Bayesian updating algorithm
The BN detailed in Figure 2 is used to perform a predictive analysis, starting from the prior distribution of the root nodes and a diagnostic analysis, entering an observation at the nodes IM MS2 , IM MS3 , RD obs and PA obs and visual observations (VI). For this purpose, the OpenBUGS software (e.g., 53 ) is employed, which is interfaced with the R statistical tool. Open-BUGS can treat both deterministic (e.g., M MS and ) and probabilistic (e.g., IM MSi , TD) variables. These latter are sampled through a Markov Chain Monte Carlo (MCMC) sampling scheme. Each chain is built with a Gibbs sampling scheme, where variables are sampled successively from the posterior distribution of previous variables: the posterior distribution of a variable is obtained from the product of the prior distribution and the likelihood function (probability of a given observation occurring given the prior distribution). The samples are then aggregated to estimate empirical statistics of the variables of interest, representing the posterior distributions. Although Bayesian inference based on sampling provides only approximate solutions (i.e., the posterior distribution is built from the samples), it has the benefit of being much more flexible than exact inference algorithms such as junction-tree inference (e.g., 54 ), since it allows modelling continuous variables using various probability distributions. Due to the approximate nature of the posterior distributions sampled by the MCMC scheme, there is no guarantee that exact distribution parameters may be obtained. However, in the present study, the following steps are taken to ensure reasonable accuracy of the results: • generation of multiple MCMC chains starting with different combinations of initial conditions to ensure that all chains converge towards the same values; • generation of a high number of samples for each chain (e.g., several tens of thousands); • definition of a 'burn-in phase', where the first part of each chain is removed from the estimation of the posterior distribution, to remove samples that have not yet converged; • thinning of the samples (i.e., only one sample in every five is considered in each chain) to reduce autocorrelation effects inherent in MCMC sampling.
Specific statistical tools in OpenBUGS are dedicated to estimating autocorrelation and require a minimum number of samples. In any case, preliminary tests are necessary to calibrate the sampling parameters carefully. The chosen sampling results from a trade-off between the required accuracy level and the computational cost. In the present application, the relatively modest size of the BN does not imply unreasonable computational times, and the convergence towards an accurate posterior distribution is checked by estimating the 'R-hat' value (e.g., 55 ) and the effective sample size (e.g., 56 ) for all variables.

Case-study description
For demonstration purposes, the structural system considered in this study consists of a two-span bridge with a continuous multi-span steel-concrete composite deck, hypothetically located in L'Aquila, Italy (latitude 42.5650N, longitude 12.6438E). The selected bridge represents a class of regular medium-span bridges commonly used in transportation networks (e.g., [57][58][59] ) (see Figure 3). The bridge superstructure, designed according to the specifications given in Eurocode 4 (ECS), 60 59 A three-dimensional finite element (FE) model of the bridge is developed in OpenSees (e.g., 61 ) following the same approach as described in Tubaldi et al., 62 that is, using linear elastic beam elements for describing the deck, and the beam element with inelastic hinges developed by Scott et al. 63 to describe the pier. Further details of the FE model and the pier properties are given in Tubaldi et al. 62 The elastic damping properties of the system are characterised by a Rayleigh damping model, assigning a 2% damping ratio at the first two vibration modes. The FE model described in this study is assumed to be deterministic and characterised by no epistemic uncertainties. Future extensions of the methodology will consider the introduction of modelling uncertainty (such as considering the approach outlined in Tubaldi et al. 64 and their effects on the results). Figure 4 shows the hysteretic response of the pier to a bi-directional ground-motion record, in terms of momentcurvature of the base and base shear-top displacement, along the two principal directions of the bridge. It can be observed that some degradation of stiffness and pinching characterise the model. This results from the constitutive model adopted to describe the concrete fibres in the plastic hinge region (Concrete 02 in OpenSees; (e.g., 61 )). However, a more sophisticated description of the hysteretic behaviour of the pier and other bridge components (see for instance 65,66 ) is out of the scope of this study.
To develop the probabilistic seismic models that describe response and damage under the MS-AS sequences, 200 MS-AS recordings from the database developed by Goda and Taylor 34 and Goda et al. 67 are considered. The database contains 703 MS events and their corresponding strongest ASs (described in terms of peak ground acceleration and magnitude). The M-R rup distribution of the data used is provided in Figure 5, along with corresponding histograms. In addition, Appendix A provides further details on the selection process of the 200 MS-AS sequences employed in this study to build the probabilistic seismic response and damage models. The RotD50 pseudo-spectral acceleration at the fundamental period of the bridge (RotD50Sa) is the selected IM for both the MS and the ASs for simplicity. More advanced IM could be used (as done in the record selection, see Appendix A) but this is outside the scope of this paper. The GMM of Lanzano et al. 68 is used to estimate the ground motion IM values at the site from the considered earthquake point sources, assuming soft soil conditions (V s = 300 m/s) and a normal fault mechanism.
As described in Appendix A, the ground motion components are scaled such that their RotD50Sa matches the target spectral acceleration level. After the ground motion components are selected and scaled, the two ground motion components are randomly rotated and then applied to the two orthogonal directions of the bridge structure. This is done to avoid any assumptions regarding the incident angle of the ground motions with respect to the bridge structure (e.g., 69 ). The PSDM described in Subsection 2.3 is built using the 200 samples of the various response parameters of interest for the performance assessment obtained under the MS event. These are the RD (EDP 1 = RD), the peak TD (EDP 2 = TD), the PAs (EDP 3 = PA), the maximum strain of concrete cover under compression (EDP 4 = ), the maximum strain of concrete cover under tension (EDP 5 = ) and the damage index D MS . The value of the coefficient considered in the calculation of the damage index is 0.8. A sensitivity analysis on the number of samples has been performed and while the analysis is not reported in the paper for brevity, the results have shown that the PSDM parameters and also the risk estimates shown in the later section do not change significantly by increasing the number of samples. Figure 6 shows the sample values of these parameters versus IM MS in the log-log plane. In the same figures, the median of the fitted PSDM is also plotted. For simplicity, the same value of IM* is used for the various parameters of interest. The value of IM* = 7.39 m/s 2 corresponds to a median value of the magnitude of the top displacement vector of 0.0106 m. According to Figure 4B On the other hand, the correlation between TD and PA is of the order of 0.7 for the first branch of the PSDM, and of 0.497 for the second branch of the PSDM. This suggests that the information on accelerations may be used to reduce uncertainty in estimating the bridge's peak TDs. It is noteworthy that this approach avoids the need to doubly integrate the measured acceleration signal when estimating displacements, which is characterised by several limitations (e.g., 70 ).
The model for damage accumulation corresponding to Equation (12) Figure 7 shows the samples of the damage index after the AS as a function of the IM of the AS and the damage index at the end of the MS , and the surface corresponding to the fitted median model. The coefficient of determination of the model, R 2 , is 0.697. These values reveal a relatively adequate model fit to the generated damage index data. To increase the accuracy of the results, a higher-order nonlinear model would be required. The regression coefficients for the models of Equations (7) and (12), the covariance matrices of the PSDM for the MS damage assessment and the corresponding correlation matrices can be found in Appendix B.
It is assumed that the bridge is equipped with one accelerometer mounted at the superstructure level above the pier. The measurement error of the accelerometer is characterised by a normal distribution with zero mean and a standard deviation of 0.002 m/s 2 . This value is based on the noise root mean square (RMS) levels of exemplary low-cost sensor specifications extracted from representative datasheets ([e.g., 71 ] for a typical low-cost MEMS accelerometer). Other sources of information are not considered in this specific study.

Rapid damage assessment for a single scenario
This subsection describes the results of the Bayesian updating of the AS risk for a single MS scenario, considering the information provided by an accelerometer mounted on the bridge deck above the pier and by a VI carried out at the end of the MS. It is assumed that seismic stations are not sufficiently close to the site to provide valuable information for hazard and risk updating. The use of information from stations close to the site is explored in the previous paper underpinning this one (e.g., 21 ). The seismic scenario corresponds to a seismic event originating from a point source with magnitude M MS 6.5, located 15 km from the site. These values are consistent with the modal values of the seismic hazard disaggregation for the region of interest. The predictive analysis is first run based on the information at the root nodes. Subsequently, multiple independent diagnostic analyses are performed by entering a piece of evidence one at a time at the nodes PA obs and by entering all the information at these nodes at the same time. These analyses are performed with OpenBugs (e.g., 53 ) using three MCMC chains generated with different combinations of initial conditions. This is to ensure that the three different starting points converge towards similar posterior distributions. Each chain contains 10,000 samples obtained by starting from 60,000 iterations, discarding the first 10,000 (burn-in), and thinning to reduce autocorrelation. Ultimately, a total of 30,000 samples are used to estimate the posterior distributions. It is noteworthy that the computation time required to perform a single Bayesian inference analysis is low (in the order of a few seconds on a standard personal computer). Figure 8 shows the empirical cumulative distribution function (ECDF) of the prior distribution of the damage at the end of the MS event D MS , and the posterior distributions given the observations of the VI ( Figure 8A) and of an accelerometer on the bridge deck ( Figure 8B). In the analyses, and have been assumed as lognormal random variables with mean values respectively equal to 0.004 and 0.001, and a coefficient of variation of 0.3. The visual inspector has been assumed to be well trained, so that concrete crushing and cracking are always correctly detected if they occurred. First of all, it is observed that the expected value of the damage index at the end of the MS is low according to the prior estimate of damage. However, the outcomes of the VI can change this distribution significantly. Two different outcomes of the VI are considered. The first one corresponds to the observation of concrete cracking and crushing, while the second one observes no cracking or crushing. In the case of the second outcome, the ECDF shifts to the left (i.e., the damage is overall less than the prior estimate), whereas in the case of cracking and crushing, it shifts significantly towards the right (damage significantly higher than expected). Two different observations of the accelerometer placed at the deck level are considered for updating the damage. In the case of a low value of recorded acceleration (PA = 2.95 m/s 2 ), the expected value of damage reduces, whereas it increases for a very high value of the recorded acceleration (PA = 7.83 m/s 2 ), as expected.
A series of 10,000 MS-AS sequences is generated following the methodology outlined in Subsection 2.2. Each sequence is described by the seismic intensities of the ASs that occur within a time window of interest following the MS, collected in the vector IM AS . Figure 9 shows the ECDF of the number of AS occurrences in an interval of 10 and 360 days following the MS. Obviously, the number of occurrences can assume only discrete values. It can be observed that the average number of occurrences of ASs of any intensity increases for increasing time. The average number of occurrences is about 1 in 10 days, and it increases to 2 in 360 days. Figure 10 shows the probability of exceedance of the IM for the case of the MS event only and for the MS-AS sequence within different time windows. It can be observed that the probability of exceedance increases significantly if ASs are considered, despite the magnitude of the AS being constrained to be less of that of the MS. This is because the AS could originate from a point source closer to the site than the source of the MS, but also due to the uncertainty inherent in the GMM regarding the ground motion attenuation from source to site. Obviously, increasing the time window of interest, the probability of exceedance also increases, but the highest relative increase is observed when only a few days after the MS are considered. This is because the rate of occurrence of the AS decreases with time since the MS occurrence.
For each of the 10,000 sample MS-AS sequences, the evolution of the damage index is evaluated by using the model corresponding to Equation (13). The various samples of the MS damage index D MS obtained by running the Bayesian network are considered as starting points for the various sequences. Obviously, different sample sets must be used, depending on whether the prior estimates or the posterior estimates of D MS following an observation are considered. The IM AS samples are then used together with Equation (13) to generate samples of the damage index vector D AS collecting the damage index values at the end of each AS. The damage at the end of the MS-AS sequence, denoted as D AS,max , is then used for evaluating the AS risk. Figure 11 shows the result of the application of the procedure to a sample sequence, considering a time window of 360 days. In particular, Figure 11A shows the values of the IM of the MS (event 1) and of the subsequent ASs (events [2][3][4][5][6][7][8]. Figure 11B shows the corresponding values of the damage index at the end of the MS and the subsequent ASs. The damage increase is observed only in correspondence to events with significant intensity, as expected. The samples of D AS,max obtained for the 10,000 considered sequences are used subsequently to estimate the probability of damage exceedance in the time window of interest. Figure 12 shows the probability of exceeding different levels of damage at the end of the MS and the end of the MS-AS sequence, 10 days after the MS. In Figure 12A, the risk related to prior estimates is compared with that related to posterior estimates that result from a VI in which no cracking or crushing was observed. In Figure 12B, the risk related to prior D MS estimates is compared with that related to starting posterior estimates that result from a measurement of a low PA (2.95 m/s 2 ) at deck level during the MS. It can be observed that the increased risk of bridge damage obtained considering the ASs hazard is low, though not negligible. Moreover, accounting for the observations from VIs or accelerometers can change the damage risk estimates significantly. The decrease in damage risk due to the observation of no cracking or crushing is more significant than the increase in damage risk due to the consideration of ASs. It is noteworthy that the obtained results are not significantly affected by the number of samples used to fit the demand and damage models. Figure 13 shows the time-dependent probability of exceeding various levels of the damage index under the MS-AS sequence, obtained starting from the prior estimate of D MS or from the posterior estimate with no observed cracking/crushing. In general, the probability of exceedance increases for an increasing number of days elapsed after the MS. The rate of increase is higher a few days after the MS and then decreases due to two different effects: the time-decaying rate of ASs and the fact that the intensity of the AS is limited. A similar trend was observed in Jalayer et al. 72 It is also interesting to observe that the relative increase of risk for the MS-AS sequence compared to the MS is more significant if the posterior estimate of D MS is considered as an initial condition. Starting from the prior estimate of D MS , the probability of D AS,max exceeding the value of 1, corresponding to significant damage, is equal to about 2.67 × 10 -2 in the case of the MS event only, and it increases to 1.22 × 10 -1 if the whole MS-AS sequence within a time window of 360 days is considered. On the other hand, starting from the posterior estimate of D MS , the probability of exceeding the damage level of 1 increases from 5.41 × 10 -6 at 0 days to 6.82 × 10 -2 after 360 days. These results could be communicated to transport agencies or bridge managers and can help make better-informed decisions concerning bridge closure in the aftermath of an MS.

CONCLUSIONS
This paper illustrated a Bayesian framework for the AS risk assessment of bridges. The main novelty aspect of the framework is that it allows to specifically exploit heterogeneous information from seismic stations, SHM sensors (accelerometers and GPS receivers) and VIs for updating knowledge of the damage state of the structure after the MS and thus make a more accurate assessment of the risk due to future ASs. The proposed framework is applied to a hypothetical bridge in Central Italy subjected to a moderately strong MS scenario. A probabilistic model is fitted to describe the joint distribution of the various parameters related to the bridge state (damage index), the monitored responses (e.g., accelerations, RDs) and the VI outcomes (concrete cracking and crushing). The correlation between the involved variables is exploited to update the damage state given the sensor readings and the VI reports.
The study results show that observations of sensors and VIs can significantly affect the decision-making process concerning bridge closure in the aftermath of an MS. For example, a bridge for which prior knowledge of the MS magnitude and source-to-site distance results in a high probability of large damage after an MS may be classified as safe following a VI that reports no cracking/crushing of bridge piers. This obviously has an impact on the risk due to future ASs. Similar results may be obtained by exploiting the information from accelerometers, if the recorded maximum absolute accelerations are low compared to the ones based on the only knowledge of the earthquake magnitude and location. It is noteworthy that the results obtained by applying the proposed framework to the case-study bridge can be strongly affected by the modelling choices. These include the description of the earthquake AS hazard, the structural model and the damage accumulation under multiple events, and the outcomes of the VIs given a possible damage state. Further improvements can be made by developing more sophisticated FE models and more consistent MS-AS ground motion selection approaches that better represent site-specific seismic hazard with a more comprehensive ground-motion database. The inclusion of more causal parameters in the analysis and more sources of information about the seismic response (e.g., information from drone-based surveys and damage detection) can also contribute to enhancing the proposed framework.
More research is needed to develop accurate models, for example, describing the relationship between amplitude and distribution of cracks and/or length of crushed concrete zone to the drifts and damage experienced by bridge piers. This will be the object of future investigations. In addition, a robust definition of the levels of damage that should trigger decisions concerning the bridge serviceability should also be the object of further study. Future works will also explore the effectiveness and benefit of alternative monitoring and inspection strategies in terms of better-informed decisions concerning bridge operability, using concepts such as expected utility theory, multi-criteria decision making and value of information.

A C K N O W L E D G E M E N T S
This paper is supported by the European Union's Horizon 2020 research and innovation programme under grant agreement No 821046, project TURNkey (Towards more Earthquake-resilient Urban Societies through a Multi-sensor-based Information System enabling Earthquake Forecasting, Early Warning and Rapid Response actions). Input and feedback from Gemma Cremen and Dr Elisa Zuccolo on the draft manuscript is greatly appreciated.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX A
The 703 mainshocks (MS) -aftershock (AS) ground-motion (GM) sequences selected from the database provided by Goda et al. 34 and Goda et al. 67 are used to compute their RotD50Sa spectra (e.g., 73 ) for 2% damped oscillator at 40 periods. The 703 MSs and AS spectra are presented in Figure A1. To select the ground motions appropriately while considering higher mode effects and period elongation effects, average RotD50Sa (denoted as RotD50Sa,avg) is utilised as the selected groundmotion intensity measure. RotD50Sa,avg is computed as the mean of RotD50Sa between 0.5 T1 and 2 T1 (e.g., 74 ) where T 1 = 0.432 s. The histograms of the RotD50Sa,avg computed for the database, are presented in Figure A2.
Firstly, MS GMs are chosen such that 10 GMs are selected for 20 levels of RotD50Sa,avg ranging from 0.1 to 2.0 g at an interval of 0.1 g with the requirement of no scaling. This leads to a total of 200 MS unscaled GMs. Figure A3A shows the histogram of RotD50Sa,avg of the selected MS unscaled GMs. It can be observed from the histogram, ∼10 GMs are selected for each level of the 20 RotD50Sa,avg with slightly lesser number of GMs available for higher RotD50Sa,avg (i.e., RotD50Sa,avg > 1.6 g). This due to the lack of recorded GM options for higher intensity levels. Finally, the 200 unscaled GMs are scaled such that 10 GMs exactly match each level of the RotD50Sa,avg to have precisely 10 GMs for each of the 20 RotD50Sa,avg levels. The histogram of the scaling factors is shown in Figure A3B, and it can be observed that they range very close to 1, thereby signifying lower scaling levels.
The natural unscaled AS GMs corresponding to the 200 unscaled MS GMs selected above are selected for the AS sequence to maintain MS-AS IM correlations. This leads to 200 unscaled AS GMs whose RotD50Sa,avg is presented in Figure A4A. Figure   the corresponding natural 200 unscaled AS GMs. It can be observed that for lower levels of MS RotD50Sa,avg levels, the RotD50Sa,avg levels of AS GMs tend to be very small, and the scatter plot tends to be very sparse for high AS RotD50Sa,avg levels. To remedy this, the selected unscaled AS GMs are scaled (with minimal scaling) such that for each MS RotD50Sa,avg level (i.e., 20 levels ranging from 0.1 to 2.0 g at an interval of 0.1 g), the corresponding ten unscaled AS GMs lead to ten levels of RotD50Sa,avg ranging from 0.1 to 1.9 g at an interval of 0.2 g. The histogram of the utilised scaling factors is shown in Figure A4B, and the RotD50Sa,avg of the selected scaled MS and AS GMs is shown in Figure A5B. To ensure that minimal scaling is used for the AS GMs, for each level of MS RotD50Sa,avg, the corresponding ten unscaled AS GMs are scaled such that each one requires minimal scaling to achieve RotD50Sa,avg between 0.1 and 1.9 g. This process leads to MS-AS sequences with RotD50Sa,avg shown in Figure A5B.

MS-AS GMs
Furthermore, the magnitudes of the MS and AS remained consistent with the statistical knowledge obtained from the history of recordings where the selected MS-AS led to mean MS magnitude = 8.15 and mean AS magnitude = 6.80. This is consistent with the literature available that mentions that the magnitude of AS is on average 1-1.2 points lower than the magnitude of MS event (e.g., 75 ). The magnitudes of the MS and AS events of selected ground motions are shown in Figure A6A and B, respectively. To further introduce randomness in the ground motions, the 'as recorded' components of the 200 MS and AS ground motions were rotated at random angles. However, it was made sure that both MS and AS GMs of the sequence were rotated to the same random angle in order to maintain the internal correlations.

APPENDIX B
The regression coefficients for the average model of Equation (7) are listed in Table B1.
The covariance matrices ∑ and ∑ , collecting the information on the variance of the error variables (in the lognormal space) and their correlation, for the two branches of the PSDM for the MS damage assessment (corresponding respectively  Table B2 presents the regression coefficients for the average model of Equation (12). The estimate of the standard deviation of is 0.603.