A generalized ground-motion model for consistent mainshock–aftershock intensity measures using successive recurrent neural networks

Several recent studies have investigated the risk posed to structures by earthquake sequences, proposing state-dependent fragility/vulnerability models for assets in damaged conditions. However, a critical component for such efforts, i.e., ground-motion record selection, has received relatively minor consideration. Specifically, utilization of “consistent” mainshock (MS)–aftershock (AS) ground motions is desirable in practical applications. Such consistency in selecting MS–AS sequences requires proper consideration of the correlations between and within the intensity measures of MS and AS ground motions. Most of the studies in this domain utilize spectral accelerations as the considered ground-motion intensity measures and rely on empirical linear correlation models between the intensity measures of MS and AS ground motions for developing, for instance, record selection approaches. This study proposes a generalized ground-motion model (GGMM) to estimate consistent 30 × 1 vectors of intensity measures for mainshocks (denoted as IMMS) and aftershocks (denoted as IMAS) using a framework of successive long-short-term-memory (LSTM) recurrent neural network (RNN). The vectors of IMMS and IMAS consist of geometric means of significant duration (D5-95,geom\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{5 - 95,geom}$$\end{document}), Arias intensity (Ia,geom\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{a,geom}$$\end{document}), cumulative absolute velocity (CAVgeom\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$CAV_{geom}$$\end{document}), peak ground velocity (PGVgeom\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$PGV_{geom}$$\end{document}), peak ground acceleration (PGAgeom\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$PGA_{geom}$$\end{document}) and RotD50 spectral acceleration (SaT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_{a} \left( T \right)$$\end{document}) at 25 periods for both MS and AS ground motions. The proposed RNN-based GGMM is trained on a carefully selected set of ~ 700 crustal and subduction recorded MS–AS sequences. The inputs to the framework include a 5 × 1 vector of source and site parameters for MS and AS recordings. The residuals of the trained LSTM-based RNN are further used to develop empirical covariance structures for IMMS and IMAS. The proposed framework is finally illustrated to select MS–AS ground motions based on IMMS and IMAS using a multi-criteria objective function. The selected MS–AS ground motion sequences are then used to perform non-linear time-history analyses of a case-study two-spanned symmetric bridge structure. The obtained engineering demand parameters are evaluated and critically discussed.


Introduction
Current structural design and analysis procedures in earthquake engineering continue to be primarily based on the hazard demands posed by mainshock (MS) ground motions. However, recent earthquake events and research studies worldwide (e.g., Stewart et al. 2018;Kam et al. 2010;Goda 2015;Shokrabadi and Burton 2018;Papadopoulos et al. 2019Papadopoulos et al. , 2020Shcherbakov 2021) have emphasized that sequences of damage-aggravating aftershock (AS) ground motions often follow moderate-to-large MS earthquakes. Thus, structures that sustain damage from a MS and are not immediately repaired are often prone to damage accumulation from subsequent AS events/ground motions. Although this is a widely recognized issue, this form of hazard interaction is commonly overlooked in current seismic hazard analysis, structural design, and risk assessment procedures. More in general, structures located in seismically active regions may frequently experience more than one earthquake event during their lifetime. Hence, estimates of structural performance and fragility models (i.e., probability of exceeding various damage states as a function of hazard intensity measures) often need to be updated soon after a structure is exposed to intense ground shaking due to a given event. This approach can better support inspection, revaluation, and repair decisions to ensure the structure's serviceability during future events, particularly during the post-event response/emergency phases (e.g., Jalayer et al. 2011;Franchin and Pinto 2009).
In recent years, several research studies have focused on quantifying the impact of sequences of ground motions (particularly MS-AS sequences) on seismic hazard analysis, structural performance, and risk assessment (Yeo and Cornell 2009;Cornell 2004;Luco and Bazzurro 2004;Shokrabadi and Burton 2018;Jalayer and Ebrahimian 2017). Many studies promote the concept of damage-dependent (or state-dependent) fragility models, i.e., fragility relationships for structures with pre-existing damage conditions (Aljawhari et al. 2020;Gentile and Galasso 2021;Raghunandan et al. 2015). Their numerical derivation often entails non-linear dynamic analyses in which the structure is subjected to backto-back MS-AS ground motions, selected through various approaches (e.g., Goda 2015;Papadopoulos et al. 2020). FEMA 352 (FEMA 2000), Cornell (2004), and Luco and Bazzurro (2004) were among the very first studies that investigated the effects of AS ground motions on buildings; they focused on inspecting post-earthquake damage by generating a demand hazard curve for damaged buildings subject to AS ground motions. These studies involved a basic assumption, i.e., independence between MS and AS ground-motion intensities.
To perform structural analyses using a proper hazard representation of MS-AS groundmotion features, sequences of ground motions that best represent the characteristics of MS-AS ground motions at the site of interest must be utilized. In particular, MS and AS ground-motion selection should ensure that the correlations between and within the intensity measures (IMs) of MS and AS ground motions are accurately characterized and considered (e.g., Papadopoulos et al. 2019). Studies such as Kohrangi et al. (2017), among others, have demonstrated the impact of alternative site-specific record-selection strategies on MS fragility functions. Due to the similarity in causal parameters such as source rupture properties, wave propagation paths, and site characteristics between the MS and subsequent AS ground motions, it can be substantiated that such site-specificity should also be accounted for in the selection of AS records. Apart from this, the selected AS records must also be consistent with the selected MS records in terms of IMs and their characteristics. In the current state of practice, such site-specificity is accounted for by considering casual parameters such as magnitude, rupture-to-site distance, and soil shear-wave velocity at the site, and the record consistency is quantified mainly in terms of spectral accelerations ( S a (T) ). While conventional Probabilistic Seismic Hazard Analysis (PSHA) relies on the disaggregation of the MS hazard to obtain the casual event parameters, recent methods of Aftershock Probabilistic Seismic Hazard Analysis (APSHA) and Omori's law (Yeo and Cornell 2009;Žalohar 2018) allow the computation of MS-consistent seismic event parameters for AS events. Several alternate methods available in the literature can be used to compute the seismicity/source parameters of AS events from the parameters of MS events. One of the most popular models to describe the seismicity of a region is the space-time Epidemic-Type Aftershock Sequence (ETAS) model (e.g., Ogata 1998;Iacoletti et al. 2022). However, due to the computational complexity of the ETAS model, particularly in terms of the location of AS events, more straightforward approaches such as the Branching Aftershock Sequence (BASS) model (e.g., Turcotte et al. 2007) are available to compute consistent source-to-site distances between MS and AS events probabilistically.
Using these advances, several researchers (Hu et al. 2018;Burton et al. 2017;Jalayer and Ebrahimian 2017;Aljawhari et al. 2020;Gentile and Galasso 2021;Goda 2015;Fayaz et al. 2019;Park et al. 2018) have proposed various solutions to perform structural analysis using appropriate sets of MS and AS ground-motion sequences. One of the most notable studies in this field was undertaken recently by Papadopoulos et al. (2020). In this research, based on the assumption of lognormality and linear correlations, a MS-AS conditional spectrum (MSAS-CS) was proposed, which models the joint distribution of AS spectral ordinates conditioned on the MS spectrum and the rupture characteristics. The developed MSAS-CS was then used to select AS records for performing Nonlinear Time-History Analysis (NLTHA) of two two-dimensional non-ductile reinforced concrete frame (RC) buildings. The selected MS-AS records were finally used to develop state-dependent fragility relationships of the building models, demonstrating the importance of considering spectrally consistent MS and AS records. However, an over-reliance on the sufficiency and efficiency of the S a (T) spectrum and linear correlations among the S a (T) ordinates can be detrimental to the accuracy of any derived conclusions. This was indicated by Fayaz et al. (2021b), who proposed an artificial neural network (ANN)-based generalized groundmotion model (GGMM) for MS records, which uses the seismic source and site parameters to output a vector of 29 × 1 ground-motion IMs. The GGMM framework was compared with the conditional spectrum approach (Lin et al. 2013), which demonstrated the superiority of ANNs in modeling cross-IM-dependencies compared to the linear correlations generally used in earthquake engineering.
This study is a step towards using current state-of-the-art deep learning-based frameworks in structural and seismic analysis. It incorporates a data-driven framework of successive recurrent neural networks (RNNs) to develop a GGMM for two 30 × 1 vectors of MS and AS IMs (denoted as IM MS and IM AS , respectively). In particular, to incorporate the higher-order dependencies among the components of IM MS and IM AS , long-short-term memory (LSTM) cells are adopted in the proposed RNN-based framework. The inputs to the proposed MS-AS GGMM include a 5 × 1 vector of seismic source and site parameters for the MS (denoted as ) and AS (denoted as ). The outputs of the framework (i.e., the IM MS and IM AS vectors) include geometric means of significant duration ( D 5−95,geom ), Arias intensity ( I a,geom ), cumulative absolute velocity ( CAV geom ), peak ground velocity ( PGV geom ), peak ground acceleration ( PGA geom ) and RotD50 spectral acceleration ( S a (T) ) at 25 periods for both MS and AS ground motions. The use of LSTM-based RNNs helps maintain the internal cross-dependencies between and within the IM MS and IM AS vectors while leading to good mean predictions. This ensures that the MS and AS are consistent in terms of their IMs. Furthermore, the residuals of the trained RNN are used to develop two 30 × 30 covariance matrices of IM MS and IM AS . The predictions of the developed GGMM are analyzed through rigorous statistical analyses. The application of the model is finally demonstrated by performing MS-AS ground-motion selection for NLTHA of a twospanned two-column symmetric bridge structure.
The paper starts by describing the ground-motion database and the proposed neural network-based models. The results of the statistical tests for the proposed models are discussed subsequently. Then, the covariance structures developed from the residuals of the RNN models are defined in the following section. Finally, the proposed GGMM is illustrated for scenario-based MS-AS ground motion selection and NLTHA of an ordinary bridge structure.

Ground-motion database
This study uses the MS-AS ground-motion sequences selected by Goda and Taylor (2012) and Goda et al. (2015) from the Next Generation Attenuation Relationships for Western US (NGA-West2) (Ancheta et al. 2004) and KKiKSK (2012) databases. The moment magnitude and rupture distance metadata of the selected 703 MS ground motions and corresponding AS ground motions from 40 events are provided in Fig. 1a and b, respectively. The dataset consists of 294 ground motions from crustal earthquake sources and 409 from subduction earthquake sources. Goda and Taylor (2012) and Goda et al. (2015) provide a detailed explanation of the rationale behind selecting these ground-motion sequences. In general, ground-motion sequences with mainshock magnitude ( M MS ) ≥ 5 and PGA geom > 0.075 g are selected, and among the recorded aftershocks, the ones with the highest magnitude are used.
In this study, no consideration is provided for the spatial correlation of ground motions (Jayaram and Baker 2010) or multiple recordings at the same site (e.g., Kotha et al. 2017).  Fig. 2a and b. In general, it can be observed from the two figures that the spectral shape of the MS and AS ground motions tend to be different. Specifically, the decay in the RotD50S a (T) values in the MS records are lower as compared to the RotD50S a (T) values in the AS records. Hence, it is essential that the ground-motion models are differentiated for MS and AS records, also accounting for the correlations between and within the MS and AS records.

Successive recurrent neural networks (RNNs)
RNNs (Williams et al. 1986) extend conventional feedforward ANNs that attempt to model time or sequence-dependent behavior. The underlying principle of modeling temporal sequences in RNNs is the use of recurrency in the networks. This is done by feeding back the output of a neural network layer at time step t to the input of the same network layer at time step t + 1. Since the ground-motion IM vectors represent the same ground motion characteristics, the values are dependent on each other. This can be viewed as a sequence of IMs where one value depends on the other values in the sequence. Therefore, RNN is suitable for training a data-driven model to predict cross-dependent vectors of IMs (IM).
However, conventional RNN structures still suffer from drawbacks, such as shortterm memory and vanishing gradients (Hochreiter and Schmidhuber 1996). LSTM units (Hochreiter and Schmidhuber 1996) are deployed in the RNN structure considered in this study to remedy this issue. LSTMs modulate the information flow using internal mechanisms of cell and forget gates. General details of LSTM-RNN structures can be obtained from Fayaz et al. (2021b). LSTM RNNs feed forward while keeping an internal memory to process the inputs sequences adaptively and maintain internal dependencies between all data points within the output vector. Furthermore, the recurrent nature of an LSTM RNN enables performing the same function for each input, copying and sending the data back to the network while simultaneously producing the output. Hence in this study, the LSTM-based RNNs are used to develop MS-AS GGMM for log-scaled intensity measures, including D 5−95,geom , I a,geom , CAV geom , PGV geom , PGA geom (also referred as S a (T = 0 s)) and S a (T ) at 25 periods (ranging from 0.05 to 5 secs) for MS and AS ground motions. Specifically, two RNN structures are trained, one for mainshocks (MS-based RNN) and the other for aftershocks (AS-based RNN), linked together, as shown in Fig. 3. The neural networks are trained with cross-validation using a randomly selected 80% of the total dataset. The MS-based RNN is trained using five seismic source and site parameters, including earthquake type ( EQ ), moment magnitude ( M MS ), closest rupture distance ( R rup,MS ), hypocentral depth ( Z hyp,MS ) and average soil shear-wave velocity for the top 30 m ( V s30 ) (collectively denoted as ). The parameter EQ is input as a one-hot vector with [1,0] for crustal and [0,1] for subduction sources.
The RNN is trained to output the corresponding IM MS vector. The prediction power of the trained RNN is analyzed by comparing the observed and median predicted values of the IM MS vector. The results of this comparison are presented in Fig. 4. It can be observed from Fig. 4 that the observed versus predicted data points lie very close to the identity line with no evident bias for any IMs. Figure 4 also provides the coefficient of determination R 2 of each IM prediction, which is observed to be greater than 0.7 in most cases. It should be noted that the framework is trained to possess good prediction power for each IM and maintain the internal cross-dependencies within IM MS . While the internal IM dependencies in IM MS are not necessarily linear, for the sake of interpretability and brevity, the empirical correlation matrices of the observed IM MS vector and predicted IM MS vector are compared. The two correlation matrices are presented in Fig. 5 for 10 out of 30 IMs in the IM MS vector. The observations suggest that, in general, the two correlation matrices are similar, and the RNN framework performs well in maintaining the internal correlations in the IM MS vector. The correlation values differ mainly within the −0.1 to 0.1 range, which is not expected to alter the results significantly. However, it should be noted that an RNN framework does not explicitly intend to preserve linear correlations within a sequence but instead tries to make sure that the internal dependencies are maintained, which can be highly non-linear. Hence a mismatch in linear correlations does not necessarily denote lower performance of the trained RNN. To ensure that the predictions of the AS-based RNN are consistent with the predictions of the MS-based RNN, the outputted IM MS of the MS-based RNN is concatenated with AS inputs of seismic source and site parameters to train the AS-based RNN. The AS-based RNN is also trained using the same five seismic source and site parameters: earthquake type ( EQ ), magnitude ( M AS ), closest rupture distance ( R rup,AS ), hypocentral depth ( Z hyp,AS ) and average soil shear-wave velocity for the top 30 m ( V s30 ) (collectively denoted as ). While the earthquake type (EQ) (which specifies the source as either crustal or subduction) cannot change between a MS and its AS events, this parameter is included as an input to the aftershock RNN to allow the neural networks to differentiate between the ASs of crustal and subduction sources. Similarly, in the case of V s30 , if records for MS and AS   Fig. 6. This comparison shows a slight drop in R 2 compared to the MS-based RNN; however, the data points lie very close to the identity line with no significant bias for most IMs. Notably, for lower levels of I a (< 0.05 m/s), the proposed RNN-based framework tends to overestimate the observed values slightly; predictions for such lower levels of I a can be further improved with more usable datasets. Figure 7 compares the observed and true correlation matrices for 10 out of the 30 IMs in the IM AS vector. Generally, the two correlation matrices are in good agreement and the RNN framework performs well in maintaining the internal Furthermore, to verify if the two successive RNNs maintain the correlations between IM MS and IM AS, the empirical correlation structures between IM MS and IM AS are developed from the true recordings and compared against the correlation structure computed from the RNN-framework predictions. The comparison is performed using the 2 hypothesis test proposed by Satorra and Bentler (2010) at a 5% significance level. The test failed to reject the null hypothesis that the two correlation structures are statistically similar. Hence, it can be concluded that the proposed RNN framework maintains the correlation structure between IM MS and IM AS in its predictions. Figure 8 presents the observed and predicted correlation structures between selected IMs of IM MS and IM AS . It should be noted that the correlation matrices shown in Fig. 8 are not necessarily symmetric as the IMs on the abscissa correspond to the AS, while the IMs listed on the ordinate correspond to the MS. It can be observed by comparing Fig. 8a and b that the two correlation matrices are generally similar. The correlation values differ only within the -0.1 to 0.1 range, which is not expected to impact the  Figure 8c further shows the contours of the observed correlations versus the correlations maintained by the RNN for S a (T) . It can be noted that the two correlation contours are significantly overlapped, showcasing the proposed RNN's power to capture internal dependencies. Figures 5, 7, and 8 are shown for the complete dataset as the results are similar between test and train sets.

Covariance structures of residuals
The two RNN structures can estimate the median vectors of IM MS and IM AS as described in Equations 1 and 2, respectively. The normally distributed residuals MS and AS (with zero mean vectors and Σ MS and Σ AS covariance matrices for MS and AS, respectively) are further used to construct an empirical joint MS-AS ground-motion covariance matrix for the 30 intensity measures as described in Eq. 3, which consists of four quadrants of covariance matrices Σ 11 , Σ 12 , Σ 21 and Σ 22 . In this, Σ 11 and Σ 22 represent the independent covariance matrices for the residuals of MS and AS ground motions (also expressed as Σ MS and Σ AS in Eqs. 1 and 2) and are expressed in Eqs. 4 and 5. Empirical correlations between the residuals of S a (T) at 25 periods from the two covariance structures, Σ 11 and Σ 22 , are presented in Fig. 9a, b, respectively. In these figures, it can be observed that, in general, correlations are sparser in the case of AS ground motions. Similarly, Σ 12 and Σ 21 represent the covariance between MS and AS residuals. In these equations, and respectively represent the correlations and standard deviations corresponding to the IM residuals denoted by the subscripts of the terms. Using these, the conditional covariance matrix for the AS residuals can be obtained using Eq. 7 (Papadopoulos et al. 2020). While Eq. 7 allows the development of covariance matrices for AS residuals that are conditional on MS residuals, in most performance-based earthquake engineering applications (e.g., for ground-motion record selection and modification), IMs are conditioned on a single value of intensity measure (i.e., IM*; usually S a (T * ) ) (Baker and Cornell 2006). In such cases, Eq. 8 can be used to implement the required correlation (Baker and Cornell 2006). In this study, as the median estimations of MS and AS IMs are internally correlated (both within and between MS and AS) due to the use of the proposed RNN framework, for the sake of brevity, Eq. 8 can be used to independently cross-correlate the IM uncertainty bands ( ± ) of MS and AS. Fig. 9 Correlations of S a (T) residuals for a mainshocks (Σ 11 ) ; and b aftershocks (Σ 22 )

An illustrative application of the proposed GGMM
This section illustrates an application of the proposed MS-AS GGMM model for consistent MS-AS record selection and NLTHA of ordinary standard bridge structures. The model is implemented to select 30 ground-motion records for eight scenarios used for NLTHA of a symmetric two-spanned two-column ordinary box-girder bridge structure, denoted as Bridge B. The bridge consists of two equal spans, each 33.6 m in length and 23 m in width, two columns of a radius of 0.84 m and a height of 6.7 m, consisting of ∼2% longitudinal reinforcement and a fundamental period ( T 1 ) of 0.83 s. The 3D finite-element model of the bridge was developed in Openseespy (Zhu et al. 2018). The model comprises seat-type abutments, a column bent, and an elastic superstructure representing the deck. For further details on the bridge structure, see (Fayaz et al. 2020a(Fayaz et al. , 2021a.

Record selection using the proposed MS-AS GGMM
For the sake of generalization and providing a simple example, eight arbitrary MS scenarios are used in this study, including a grid-based combination of magnitudes of 6, 6.75, (4) .25, and 8, with rupture distances of 10 km and 20 km (i.e., 4 M MS values × 2 R rup,MS values = 8 scenarios) to account for different earthquake types and ground-shaking intensities. The events with M MS of 6 and 6.75 are simulated as crustal type with EQ = [1, 0] , while events with M MS of 7.25 and 8 are simulated as subduction events with EQ = [0, 1] . The magnitude and location of the MS event can be used to obtain the corresponding (and hazard consistent) AS scenarios. As discussed above, there are several alternative methods available in the literature that can be used to compute the parameters of AS using the parameters of MS.
This study primarily utilizes the general procedure adopted by Papadopoulos et al. (2019). For each of the eight scenarios, one AS scenario with a magnitude M AS is generated using Eq. 9 (Ogata 1998), and for simplicity, the rupture distance is kept the same for the MS and AS (i.e. R rup,MS = R rup,AS ). While in reality the R rup,MS and R rup,AS can be different from each other, this assumption is in line with a previous research study by Goda et al. (2012) and is not expected to alter this study's results/discussion significantly. In Eq. 9, u m is a uniformly distributed random variable over the range [0,1]; b is the Gutenberg-Richter (GR) constant parameter; and M min is the minimum considered magnitude (set to 5 in this study). For illustrative purposes, a value of b = 0.99 is used from the ETAS calibration of the Uniform California Earthquake Rupture Forecast, Version 3 (UCERF3; Field et al. 2014) for California. Table 1 presents the magnitudes and distances of the MS and AS events for the eight considered hazard levels. For the events simulated as a crustal events (i.e., EQ = [1,0]), Z hyp is randomly simulated to be between 0 and 20 km, while for events simulated as subduction events (i.e., EQ = [0,1]), Z hyp is randomly simulated to be between 20 and 50 km. For all the eight scenarios, V s30 =360 m/s is used for illustrative purposes. The source parameters of the simulated scenarios are shown in Table 1. It should be noted that this example is used only as a demonstration to utilize the proposed MS-AS GGMM for record selection and NLTHA. The model can be easily used for any other sophisticated type of analysis.
The obtained source parameters EQ, M MS , M AS , R rup,MS , R rup,AS , Z hyp,MS , Z hyp,AS and V s30 are appropriately used as and for inputting into the GGMM. The outputted and are used as the target spectra for record selection. Figure 10 showcases the constructed and for two simulated scenarios. The mean and vectors are used to select 30 MS and AS ground motions for each scenario, as described below.
To avoid any significant bias, MS and AS ground motions are selected using the respective MS and AS datasets (Burton et al. 2017;Kohrangi et al. 2017). Minimal scaling is used with scale factors between 0.5 and 2.0 (Fayaz et al. 2021a). It should be noted here that the  scaling process alters all the elements of and vectors (i.e., PGV geom , PGA geom and S a (T) are scaled linearly and I a,geom is scaled by the square of the scaling factor) of the ground motions except for the duration D 5−95,geom . Hence, to select the ground motions that match the IM vectors, firstly, the ground motions characterized by D 5−95,geom values between 80 and 120% of the target D 5−95,geom for the given scenario are used as the groundmotion selection pool from the respective datasets. Subsequently, the pooled ground motions are scaled from 0.5 to 2.0 (i.e., 0.25 to 4 for I a,geom ) at an interval of 0.1 (hence 16 scaling factors), and the resultant 29 × 16 IM vectors (without D 5−95,geom ) are obtained and compared against the target mean GGMM IM vector of the given scenario. Due to the different scales of the IM values, the comparisons are performed using a weighted average accuracy, where accuracies for I a,geom , CAV geom , and PGV geom are computed using the normalized accuracy given in Eq. 10, and the accuracy of the S a (T) spectrum is computed using the Index of Agreement (Willmott et al. 1985) in Eq. 11. Both Eqs. 10 and 11 are bounded between 0 and 1, where 0 is the worst match and 1 is the best match. Then the accuracy measures are combined using a total weighted average accuracy ( Acc avg ), which is computed as per Eq. 12. For illustrative purposes, the weights are chosen as w 1 = 0.2; w 2 = 0.2; w 3 = 0.2 and w 4 = 0.4; however, users can freely choose weights based on their own criteria and their specific applications. Also, in this study, only the S a (T) values within the period range of 0.5 T 1 (i.e., 0.41 s) and 2 T 1 (i.e. 1.66 s) (Fayaz et al. 2021a) are used in Eq. 11 for the bridge record selection. For each scenario, the ground motions and scaling and : a using scenario #7 for S a (T) ; b using scenario #7 for other IMs; c using scenario #3 for S a (T) ; and d using scenario #3 for other IMs factor combination led to the highest 30 Acc avg are selected for the bridge assessment. In addition, the selected ground motions only contain one scaled version. The process is performed for all eight scenarios for both MS and AS. Figure 11 presents the mean and of the selected 30 MS and AS recorded ground motions for scenario #3 with M MS = 7.75 and M AS = 5.07. While in this exercise, only one is obtained for each MS ( ), for more sophisticated analysis, a multitude of can be easily obtained using the GGMM by keeping the constant and appropriately changing and then selecting ground-motion records. (10)

NLTHA results and discussion
The selected MS and AS ground motions are concatenated to develop the sequence of ground motions. To allow the bridge structures to return to steady-state response after the MS, the ground motion time-history is padded with zeros equivalent to 25 times the bridge first mode period, i.e., 250 × 0.8 secs = 200 secs (Fayaz et al. 2019). The zero-padded MS ground motion is then concatenated with AS ground motion and used to conduct NLTHA of the bridge B structure for each scenario. Since the behavior and performance of bridge structures is not similar in all axes, NLTHA is conducted by rotating the two orthogonal components of the ground motions from 0 to 180 degrees at intervals of 15 degrees (13 intercept angles). For each intercept angle of the bi-directional ground motion loading, its corresponding maximum column drift ratio (CDR) value (Fayaz et al. 2020a) through the time history is obtained. Three types of maximum CDRs are recorded from the time history of the MS + AS sequence. As shown in Fig. 12a, the three CDRs include (i) maximum CDR recorded during the MS (denoted as CDR MS ); (ii) maximum CDR recorded during the AS (denoted as CDR AS ); and iii) maximum CDR recorded during the MS + AS sequence (denoted as CDR seq ). The maximums of the 13 maximum CDR MS , CDR AS , and CDR seq from 13 rotation angles are termed as Rot100CDR MS , Rot100CDR AS , and Rot100CDR seq , respectively. Figure 12b presents the Rot100CDR seq for the eight used scenarios. The cases with Rot100CDR seq > 0.1 are classified as collapses in the analysis. The bridge structure tends to yield around ~ 1% CDR. As shown in Fig. 12b, the bridge structure remains in the elastic zone for most ground motion sequences for the eight scenarios. However, there are multiple cases where the bridge goes to non-linearity and large drifts.
Since peak engineering demand parameters do not necessarily increase monotonically with the ground-motion duration/sequence (e.g., Fayaz et al. 2020b;Gentile and Galasso 2021), it is possible that CDR MS maybe equal to CDR seq and an AS does not lead to an increase in the peak demands (unless the finite-element model appropriately captures damage accumulation and structural degradation). Figure 13 shows the comparison of Rot100CDR MS and Rot100CDR AS for the eight considered scenarios. The (12) Acc avg = w 1 Acc I a,geom + w 2 Acc CAV geom + w 3 Acc PGV geom + w 4 IA S a w 1 + w 2 + w 3 + w 4 Fig. 12 a CDRs recorded during a sample MS + AS sequence; b Rot100CDR seq for the eight scenarios figure shows the response values through scatter plots; furthermore, the kernel density estimates (KDEs) of the two types of responses are statistically tested for similarity through the Kolmogorov-Smirnov (KS) hypothesis test (Massey 1951). The p-values < 0.05 indicate that the test rejects the null hypothesis (i.e., the two samples are from the same distribution and are statistically similar). As can be observed from Fig. 13, the p-values for all eight scenarios are less than 0.05, thereby suggesting that the Rot100CDR MS and Rot100CDR AS are not similar. Moreover, by observing the scatter plot and histograms, it can be concluded that the Rot100CDR AS values, on average, are higher than Rot100CDR MS values. This means that, in most cases, the bridge model is appropriately able to capture the accumulation of damage in the structure, which leads to an increase in the drift response. It is worth stressing that this demonstration is not an intended novelty of the study but showcases how to efficiently utilize the proposed MS-AS GGMM to select and scale ground motions based on a vector of IMs and conduct proper seismic demand analysis. It is clear that consistent record selection of MS + AS sequences is essential to perform appropriate seismic demand analysis. It can lead to statistically higher demand levels that should be adequately accounted for at the design and risk assessment stages. Methods incorporating various types of dependencies between source and site, ground-motion IMs, and other characteristics of MS-AS ground motions, like the one proposed in this study, are much needed for accurate analysis and design of structures and risk assessment. Due to the preeminence of neural network-based frameworks and their capability to interpolate and extrapolate complex data (Fayaz and Galasso 2022), such methods can be used to provide efficient tools for engineers, risk modelers, and other end users. Due to the scarcity of usable earthquake recordings (especially AS records), the framework is currently trained with only ~ 700 relevant ground motions. However, the data-driven nature of the framework means it can be easily re-trained using more extensive databases for developing more generalized results. With more Fig. 13 Comparison of Rot100CDR MS and Rot100CDR AS for the eight scenarios comprehensive data, the proposed framework can also be trained on a regional basis, which may assist end-users in performing more accurate aftershock probabilistic seismic hazard and risk analyses and also provide better models for simulating synthetic ground-motion sequences.

Conclusions
This study presents a data-driven Generalized Ground Motion Model (GGMM) to predict two consistent vectors of ground motion intensity measures for mainshock and aftershock (denoted as IM MS and IM AS ) with 30 components. For this IM vector estimation, the GGMM incorporates two LSTM-based RNN structures (one for mainshock and the other for aftershock) that use two 5 × 1 vectors of inputs of earthquake source and site parameters, including earthquake type, moment magnitude, closest rupture distance, hypocentral depth, and average soil shear-wave velocity which collectively describe the physics of the rupture and characteristics of the site. The suggested IM MS and IM AS include geomeans of I a , CAV, PGV, PGA, and D 5-95 , which comprise five out of 30 components of the suggested IM MS and IM AS obtained from pairs of two horizontal components of ground motion records. The other 25 components of IM MS and IM AS include 5% damped RotD50S a at 25 periods (ranging from 0.05 to 5 secs). The suggested GGMM for estimation of IM MS and IM AS is a step in the direction that incorporates both time and frequency domain IMs of mainshock and aftershock ground motion along with their cross and internal dependencies in a single stand-alone model. The residuals obtained from the LSTM-based RNN models for the 30 IMs are further utilized to develop three covariance matrices, i.e., covariances for within-IM MS , within-IM AS , and between IM MS and IM AS . Based on the results, it is observed that the proposed GGMM demonstrates acceptable performance in predicting the vectors of intensity measures for mainshock and aftershock sequences while maintaining their internal correlations.
Finally, the use of the proposed GGMM is illustrated for MS-AS record selection for time-history analyses of an ordinary bridge structure. The bridge structure is represented using an advanced non-linear finite-element model. The proposed GGMM is used to select 30 MS-AS ground motions with minimal scaling (i.e., scale factors between 0.5 and 2) for eight scenarios. A multi-criteria objective function is used to select the MS-AS ground motions that match the target IM vectors for MS and AS. The presented example is intended to serve as a reference for utilizing the proposed GGMM for structural design and analysis, risk assessment, and decision-making.
In conclusion, the GGMM proposed in this study offers a robust data-driven tool for obtaining consistent MS and AS ground motion intensity measures and can be used for several purposes. These include pre-event structural and geotechnical design and analysis (e.g., for ground motion selection using multi-objective IM criteria), post-event riskand reliability-based decision-making, validation of artificial/simulated ground motion sequences (by checking if the simulated ground motions comply with a set of possible IM vector corresponding to causal parameters), mainshock-consistent aftershock probabilistic seismic hazard analysis, etc. Furthermore, the proposed framework can be easily re-trained with other ground motions records (if available) or extended to a larger vector of IMs that can include any other IMs of interest.