Different parameterizations, related to mixed-effects structure, of ground-motion mod- els for the duration of strong shaking are revisited. A new parameterization is proposed which can better account for the separation of source and path duration in additive ground-motion duration models. In addition, different distributions for the likelihood of ground-motion duration given the predictor variables are investigated. Traditionally, du- ration has been modeled as lognormally distributed, which makes sense for a multiplicative model, but less so for an additive model. Models using a lognormal and Gamma likelihood are compared, using the M9 simulations (Frankel et al., 2018) as an underlying data set. In general, differences between models in terms of predictions are small, but the models using a Gamma distribution perform slightly better than the ones employing a lognormal distribution. Furthermore, the new parameterization outperforms traditional models. All models are estimated via Bayesian inference, accounting for epistemic uncertainty through the posterior distribution of the parameters. 1 Introduction Typically, seismic hazard analysis is calculated form ground-motion parameters that character- ize the amplitude of the shaking. A complete description of ground motion includes a measure of duration. Ground-motion duration can influence structural damage (e.g. Chandramohan et al., 2016), liquefaction (e.g. Green et al., 2020), and slope-stability (e.g. Wang et al., 2021). Sev- eral different duration measures have been proposed, see Bommer and Mart´ ınez-Pereira (1999) for a comprehensive overview. Different duration measures can be more useful for different applications. Similar to ground-motion models (GMMs) for peak parameters or response spectral ordi- nates (pseudo-spectral acceleration, PSA), one can derive empirical GMMs for ground-motion duration measures (e.g. Bommer et al., 2009; Kempton and Stewart, 2006; Afshari and Stew- art, 2016; Jaimes and Garc´ ıa-Soto, 2021; Bahrampouri et al., 2021; Anbazhagan et al., 2017; Yaghmaei-Sabegh et al., 2022; Du and Wang, 2017). PSA GMMs model the target variable as a multiplication of different terms corresponding of source, path, and site effects. By contrast, duration GMMs come in two distinct fashions: multiplicative, where total duration is modeled as source duration times path duration (e.g. Bommer et al., 2009; Jaimes and Garc´ ıa-Soto, 2021; Bora et al., 2014; Yaghmaei-Sabegh and Hassani, 2020), and additive, where total dura- tion equals source duration pluspath duration (e.g. Abrahamson and Silva, 1996; Bahrampouri ∗kuehn@ucla.edu †mwalling@geoengineers.com 1 40.042.545.047.5 −127.5−125.0−122.5−120.0 Longitude ( °)Latitude ( °)Figure 1: Location of M9 stations used in the analysis. Red line is a contour line of RRUP= 10km. et al., 2021; Kempton and Stewart, 2006; Afshari and Stewart, 2016). An additive model as- sumes that the dominant contributor to path duration is the time difference of arrival times of compression and shear waves (Bahrampouri et al., 2021), in which case the path duration should be magnitude-independent. Multiplicative duration models have the same structure as PSA GMMs, and can be fit in the same manner, to account for between-event and site-to-site residuals via random effects models (Stafford, 2014). Additive models, however, require more care to properly account for event terms associated with event-specific source durations (Walling et al., 2018). Here, we discuss the implications of different parameterizations of additive duration models, and propose a model structure that better accounts for the separation of source and path durations than in previous models. In addition, we investigate the Gamma distribution as a way to model the likelihood of ground-motion durations. Duration models typically assume a lognormal distribution, which is very reasonable for multiplicative models, but its application to an additive model of positive terms is less natural. We fit models with both a lognormal and Gamma likelihood, and compare the estimated models in terms of their predictions. We estimate the different models based on data from the M9 simulations (Frankel et al., 2018). This provides a great data set to illustrate the different models. Compared to real data sets with observed data, the data distribution is more homogeneous. In particular, there is abundant data at short distances, which is important for duration models, since the source duration represents the model prediction at zero distance. The paper is organized as follows: we very briefly describe the data, then lay the foundation of the duration models in Section Intuition of Duration Model. Details about the M9 models are found in Section Duration Models for M9 Simulations, followed by results and a discussion. 2 Data We illustrate the different duration models on simulated ground motions from the M9 project (Frankel et al., 2018). The M9 simulations comprise simulations of 30 different magnitude 9 events on the Cascadia subduction zone. The nice thing about simulated data is that the distance range is “complete”, i.e. we do not have to worry about missing data in particular at short distances. We use a subset of the M9 simulations (mainly to reduce the number of data). Figure 1 shows a map of the simulated stations that we use. Selection is also guided by the desire to have only stations in one direction of the rupture (i.e. only towards the coast). In total, we use 15,810 observations from 30 simulations, i.e. 527 observations per simulation. 2 050100150 0 100 200 RRUP (km)D5−75 110100 1 10 100 RRUP (km)D5−75Figure 2: Used D5075duration values of the M9 simulations against rupture distance RRUP, color coded by simulation. The black line is a loess fit (Cleveland, 1979) to the data. Each simulation is treated as an event, with an event term. The predictor variable is the closest distance to the rupture RRUP, also called rupture distance. As the duration measure, we use the significant duration (Trifunac and Brady, 1975), which which is based on the normalized arias intensity I(t). Here, we use D5−75, which is commonly used (e.g. Bahrampouri et al., 2021; Kempton and Stewart, 2006; Afshari and Stewart, 2016; Jaimes and Garc´ ıa-Soto, 2021). 3 Intuition of Duration Model Generally, there are two approaches to model duration. The first is to model duration as a multiplicative model, similar to PSA GMMs (e.g. Jaimes and Garc´ ıa-Soto, 2021; Anbazhagan et al., 2017; Meimandi-Parizi et al., 2020; Bora et al., 2015, 2019). Such an approach has the advantage that the model can be fit in the same way as a response spectrum GMM (using mixed- effects regression to partition into between-event and within-event residuals, corresponding to source and path duration). A multiplicative duration model has the basic form lnD=f(⃗ r;M, R, . . . ) +δB+δW (1) where f(⃗ r;M, R) is a function with coefficients ⃗ r, dependent on magnitude, distance, and possibly other parameters such as time-averaged shear wave velocity in the upper thirty meters VS30.δBandδWare the event term and within-event residual, respectively. These are typically (universally) assumed to be normally distributed with standard deviations σSandσP. Thus, the overall duration is assumed to be lognormally distributed in a multiplicative model. A different approach to model duration is an additive model, where the total duration is modeled as the sum of a source and path duration, Dtotal=Dsource +Dpath(e.g. Afshari and Stewart, 2016; Kempton and Stewart, 2006; Boore and Thompson, 2014; Abrahamson and Silva, 1996; Bahrampouri et al., 2021; Boore, 2003; Atkinson and Boore, 1995). Empirical duration models based on regression analysis using such a model sometimes try to account for between-event terms (ostensibly corresponding to a source term) with the following modeling approach (Afshari and Stewart, 2016; Kempton and Stewart, 2006; Bahrampouri et al., 2021; Abrahamson and Silva, 1996) lnD= ln [ fsource(⃗ c;M) +fpath(⃗ r;R)] +δB+δW (2) where fsource(⃗ c;M) and fpath(⃗ r;R) are functions for the source and path duration, respectively, andδBandδWare the between-event and within-even residual. The model of Equation 2 3 has the advantage that it can be fit using a traditional mixed-effects regression (albeit with a nonlinear functional form). However, Equation (2) implies D=fsource(⃗ c;M)eδBeδW+fpath(⃗ r;R)eδBeδW(3) so the “full” source and path duration are both affected by the event term δBand within-event residual δW. In particular, the event term δB, which is typically interpreted as a term that determines the variability of the source duration, has the effect of changing the path duration (and thus the distance scaling). This means that δBcan be understood as a random effect that models an individual distance scaling for each event. Since in a model like Equation (2) the residuals δWandδBaffect both source and path durations, we call this model a “mixed model”. One way to separate the source and path duration in an additive model would be to write it as D=fsource(⃗ c;M)eδB+fpath(⃗ r;R)eδW(4) in which the two terms are completely decoupled. Assuming δBandδWare normally dis- tributed, this means that both source and path duration are lognormally distributed. The sum of two logarithmic distributions is not a logarithmic distribution, and as a consequence, the total duration Dis not lognormally distributed. Such a model was used in Walling et al. (2018). A (rather severe) constraint in Equation (4) is that all terms are positive. It then follows that the source duration, fsource(⃗ c;M), must be smaller than the smallest observed duration, D. This is a severe constraint on the source duration, and ignores possible sources of random variation in observed duration, in particular at shorter distances (see for example Figure 4(b) of Boore and Thompson (2014), which shows negative path durations after correcting total durations using a reasonable, physics based source duration model). To alleviate the constraint that the source duration must be smaller than the observed duration in an additive model, we propose to partially decouple the two terms, and write the model as D= fsource(⃗ c;M)eδB+fpath(⃗ r;R) eδW(5) In this model, the path duration is unaffected by the between-event residual δB(which is thus a proper measure for the variability of the source duration), while the structure of the model allows for the observed duration Dto be smaller than the event-specific source duration, which is given by Dsource =fsource(⃗ c;M)eδB, in the case that δWis negative. In this model, the total duration Dis lognormally distributed (assuming δWis normally distributed). We can write the model of Equation (5) as a hierarchical model Dsource∼LN(lnfsource(⃗ c;M), σS) (6) med path=Dsource +fpath(⃗ c;R) (7) D∼LN(ln med path, σP) (8) One can understand these Equations as follows: the source duration is lognormally distributed, with a median which is a function of magnitude, and standard deviation σS. Given a value of the source duration for a specific event, the median path duration is the sum of the source duration and a function that depends on distance. The actual observed duration is then lognormally distributed with the calculated median and standard deviation σP. One should note some differences of the model of Equation (5) to traditional PSA GMMs as well as the mixed model of Equation (2). The model cannot be fit using traditional mixed effects, and the total variance is not the sum of the individual variances (i.e. σ2 total̸=σ2 S+σ2 P). 4 3.1 Why Lognormal? In general, duration models assume that Dis lognormally distributed (we are not aware of any model that does not make this assumption). One reason for this assumption is ease of computation, and the advantage that one can use methods that have been well established for PSA GMMs. However, for response spectral GMMs the lognormal distribution makes physical sense, as the spectrum is a multiplication of several terms (Boore, 2003). For duration models, the situation is not as clear, in particular if one wants to use an additive model (for a multiplicative model as in Equation (1) the lognormal distribution makes sense). It is easy to reformulate Equations (6) to (8) using a different distribution, as long as it is defined for a positive random variable. Instead of the lognormal distribution, we also model duration with a Gamma distribution. In the following paragraphs, we briefly describe the Gamma distribution and how it can be used in a duration model context. If the random variate Xis distributed according to a Gamma distribution, X∼Gamma( α, β) (9) then the corresponding probability density function is pd f(x) =βα Γ(α)xα−1e−βx(10) where Γ( ·) is the Gamma function, αis called the shape parameter, and βis the rate parameter. The mean of the Gamma distribution is calculated as E[X] =µ=α β(11) If we assume that duration is distributed according to a Gamma distribution, then we model the mean as a function of magnitude (for source duration) and distance (for path duration). Hence, we can write µsource =fsource(⃗ c;M) (12) We can then rewrite the rate parameter βof the Gamma distribution as β=α µ(cf. Equa- tion (11)), and thus formulate the hierarchical model with Gamma distributions as Dsource∼G(αS,αS µsource) (13) µpath=Dsource +fpath(⃗ r;R) (14) D∼G(αP,αP µpath) (15) analogous to the lognormal model in Equations (6) to (8). It would also be possible to model the source and path duration with a different distribution (e.g. a lognormal distribution for the source and a Gamma distribution for the path duration). The Gamma model has the same number of parameters to fit as the lognormal model: co- efficients for the source and path duration, event-specific source durations, and a distributional parameter ( σfor lognormal, αfor Gamma) for both the source and path duration. One dif- ference is that we model the mean as a function of magnitude and distance for the Gamma model, while we model the median in a lognormal model. We do not think that there is an a-priori a preference for either the lognormal or Gamma model for an additive duration GMM. In the literature, duration has been universally modeled with a lognormal distribution (e.g. Kempton and Stewart, 2006; Afshari and Stewart, 2016; Bahrampouri et al., 2021; Jaimes and Garc´ ıa-Soto, 2021), but one could argue that this is due to familiarity (based on PSA GMMs), as well as ease of computation. In later sections, we compare duration models based on the Gamma model and the lognormal distribution. 5 3.2 Source Duration Model The model of the source duration is often based on the inverse of the corner frequency of the source spectrum (e.g. Abrahamson and Silva, 1996; Kempton and Stewart, 2006), following Brune (1970, 1971) Dsource =1 fc(16) Following Abrahamson and Silva (1996), the corner frequency can be written as fc= 4.9 106β∆σ M01/3 (17) where M0is the seismic moment ( M0= 101.5M+16.05),βis the shear-wave velocity at the source ( from now on we assume β= 3.2km/s), and ∆ σis the stress parameter. Thus, this give a model for the source duration as Dsource =∆σ(M) 101.5M+16.05−1/3 / 4.9 106β (18) Magnitude scaling in the source duration is introduced by the moment M0, as well as a po- tentially magnitude dependent stress parameter. Assuming an exponential dependence of the stress parameter on magnitude (∆ σ= exp( a1+a2M)), the scaling can be combined (Bahram- pouri et al., 2021; Bommer et al., 2009) into a model of the form Dsource = exp( c1+c2M) (19) A constant (magnitude-independent) value of the stress parameter ∆ σimplies a value of the slope of ln Dsource with magnitude of c2= 1.151 (Bommer et al., 2009), while a slope of a2= 3.45 for the magnitude dependence of ∆ σleads to a constant (magnitude-independent) source duration (Afshari and Stewart, 2016). We want to stress that Equation (18) is merely a crutch to connect source duration to magnitude, and to inform the form of the magnitude scaling for the source duration. While Equation (18) is developed(?) for circular ruptures, its use in empirical duration models does not mean that the events used to estimate the empirical model are assumed to have a circular rupture. Similarly, the stress parameter ∆ σis an empirical parameter that makes Equation (18) work. It should also be noted that the use of a double-corner frequency model, as used in Atkinson and Silva (2000) to model a finite fault comprised of Brune subfaults, leads to the same model structure. In this case, the source duration is modeled as (Boore and Thompson, 2014) Dsource = 0.5/fa (20) where fais the lower of two corner frequencies in the double-corner frequency model of (Atkinson and Silva, 2000), which is used to model a finite fault comprised of Brune subfaults. An empirical model for fabased on Californian events is lnfa= 2.181−0.496M (21) which leads to the same magnitude dependence as Equation (19) but with a different slope. In an empirical duration model based on events of different magnitudes the slope would be a parameter to be estimated. As stated before, in an empirical duration model one wants to account for between-even source duration variability, i.e. the fact that an event can have a longer or shorter source 6 duration than what the model implies. This can be modeled by including a random effect for each event. The distribution of the random effects allows one to include source duration variability in a forward prediction. There are different possibilities how to include an event- specific random effect in the source duration model based on the corner frequency. One can model an event-specific stress parameter ∆ σe, or an event-specific source duration Dsource,e , where subscript eindicates that the parameter corresponds to event e. If the stress parameter is modeled as a lognormal distribution, then the source duration is also a lognormal distribution, i.e. from ∆σ∼LN(µln ∆σ, σln ∆σ) (22) it follows that Dsource∼LN(µlnDS, σlnDS) (23) where µlnDSis calculated according to Equation (18) (using exp µln ∆σas the stress parameter), andσlnDS=σln ∆σ 3. This is not the case for other distributions, i.e. a Gamma distributed stress parameter does not imply a source duration that is Gamma distributed or vice versa. In general, we think that the fact that the two parameterizations of the model (using stress parameter or source duration directly) are equivalent for the lognormal distribution makes this model preferable. It is difficult to test this part of the model, as the source durations are not directly observable. The M9 simulations (obviously) have a constant magnitude M= 9, so we cannot infer a magnitude scaling. However, we can estimate a model using the stress-parameter parameteri- zation, which induces a scaling of the source duration with magnitude. For ease and stability of computation, we use the following model for the source durations of the M9 simulations lnDsource∼N(µlnDS, σlnDS) (24) i.e. with a constant median logrithmic source duration µlnDS. One can calculate a median stress parameter from µlnDSusing Equation (18). 3.3 Path Duration Model The model for the path duration needs to incorporate some sort of distance scaling that en- sures that duration increases with distance. This is often modeled as a combination of linear segments (Kempton and Stewart, 2006; Atkinson and Boore, 1995; Boore and Thompson, 2014; Bahrampouri et al., 2021; Abrahamson and Silva, 1996; Afshari and Stewart, 2016) in additive models. Based on Figure 2, we use a bilinear model for the distance scaling. The model is fpath(⃗ r;RRUP) =r1RRUP+ (r2−r1)δln 1 + expRRUP−Rb δ (25) which is a bilinear function with a smooth transition at R=Rb, where r1andr2are the slope below and above the breakpoint Rb, respectively, and δcontrols the smoothness of the transition. We fix the value of δ= 1, the other parameters are estimated during model fitting. The model of Equation (25) is differentiable with respect to all parameters, which helps during model fitting. One note on computation: for large distances, the term expRRUP−Rb δ can be numerically unstable and lead to overflow. Hence, one should use a logSumexp function, as implemented e.g. in the package matrixStats (Bengtsson, 2022) for the computer environment R(R Core Team, 2021). 7 3.4 Model Predictions Making prediction with the partially decoupled model is more complicated than for the mixed or a multiplicative model. In the mixed model, where the event term “sits outside the logarithm” in Equation (2), the median prediction can just be calculated as med= ln [ fsource(⃗ c;M) +fpath(⃗ r;R)], and the standard deviation of the corresponding normal distribution (in log-space) is calculated asσ=p σ2 S+σ2 P. By contrast, in the partially decoupled model the (event-specific) source duration is part of the prediction for the full duration, and thus there are different possibilities for calculating predictions of a new event: one could use the mean of the source duration distribution, the median, the mode, or some other point estimate. The mean of the total duration can be calculated as µ(D) =Z∞ 0Z∞ 0x p(DS)p(x|Dp)dxdD S (26) This is the same value as when using the mean of the source duration in the prediction µ(D) =µ(DS) +fpath(⃗ r;RRUP) (27) where µ(DS) = exp µlnDS+ 0.5σ2 lnDS for the M9 simulation model (as the source duration model is a lognormal distribution). This is not true for the median, i.e. med( D)̸= med( DS) + fpath(⃗ r;RRUP). If the full distribution is desired, one can use Monte Carlo simulation, i.e. one can sample from the source duration distribution, calculate desired quantities (such as mean or fractiles) of the total duration distribution given the sampled source durations, and then average over the samples. However, this approach is more computationally involved, and results may change from run to run, or depend on the number of samples. 4 Duration Models for M9 Simulations We estimate different duration models on the durations from the M9 simulations. We then compare the models in terms of their scaling, their predictive distributions, and their predictive performance. The estimated models are 1. A mixed model (cf. Equation (2)). 2. A partially decoupled model, with the source and path duration modeled as a lognormal distribution. 3. A partially decoupled model, with the source duration modeled as a lognormal distribution and the path duration modeled as a Gamma distribution. The mixed model has the form lnD= ln [ c1+fpath(⃗ r;RRUP)] +δB+δW (28) where δBandδWare normally distributed with mean zero and standard deviations σSand σP, respectively, and c1is the source duration function, which is modeled as a constant since all vents have the same magnitude. As seen in the Figure 2, there is more variability in the logarithmic duration values at short distances. Such a behavior cannot be modeled with a constant standard deviation σPor shape 8 parameter αP. Hence, we model these parameters as distance dependent, with the following functional forms σP(RRUP) =s1+s21 1 + exp( −s3(RRUP−Rs))(29) αP(RRUP) =a1+a2RRUP (30) The parameters s1,s2,s3andRs, ora1anda2, respectively, are determined during model fitting, and are constrained such that the overall standard deviaton or shape parameter is positive. 4.1 Model Implementation The models are fitted via Bayesian inference (e.g. Spiegelhalter and Rice, 2009; Gelman et al., 2013) using the program Stan (Carpenter et al., 2017; Stan Development Team, 2022). Stan implements Markov Chan Mote Carlo (MCMC) sampling (Neal, 1993) in the form of Hamil- tonian Monte Carlo (Neal, 2011; Betancourt, 2017b,a). Thus, the results of the model fitting process are samples from the posterior distribution, which can be used to assess the epistemic uncertainty associated with the models. We run the regressions through the R-interface cmdstanr (Gabry and ˇCeˇ snovar, 2021). We run four chains for each model, with 200 sampling iterations. Convergence of the chains is evaluated via the ˆR(R-hat) statistic (Vehtari et al., 2020). Comparisons between the mod- els is done via the PSIS-LOO information criterion (Vehtari et al., 2017, 2021), which is an approximation to leave-one-out cross-validation. We also assess the posterior predictive distri- bution of the models, which is a measure of the overall uncertainty/variability of the model, and comprises both epistemic uncertainty of the parameters and aleatry variability due to the lognormal or Gamma distribution. The posterior predictive distribution for an observation ˜ y is defined as (Gelman et al., 2013) p(˜y|y) =Z p(˜y|θ)p(θ|y)dθ (31) so it is calculated by integrating out the uncertainty of the estimated parameters. The parameters to be estimated are the coefficients for the path duration model ( r1,r2, Rb), the log median of the source duration ( µlnDS) for the partially decoupled model or c1for the mixed model, and the shape/scale parameters σS,σPandαP. In addition, the individual source durations DS(for the decoupled models) or the event terms δBfor the mixed model are estimated as well. The prior distributions for the source durations DSis given by Equation (24), while it is a normal distribution with mean zero and standard deviation σSforδB. For the 9 other parameters, the prior distributions are as follows: r1, r2∼N(0,0.2) RB∼N(60,40) c1∼N(25,5) σS∼N(0,1) µlnDS∼N(3.2,2) σlnDS∼N(0,2) a1∼N(0,5) a2∼N(0,1) s1, s2∼N(0,1) s3∼N(1,2) Rs∼N(60,40) 5 Results Figure 3 shows the posterior density, assessed from the posterior samples, of the parameters associated with the path duration as well as the average source duration. In general, the path duration coefficients are very similar between the three models, with slightly smaller values of r1for the Gamma model, and slightly wider uncertainty of r2for the mixed model. The average source duration is captured by parameter c1(for the mixed model) and exp( µlnDS) for the partially decoupled models (Gamma and lognormal). It models the median of the source duration distribution, and is the model prediction for RRUP= 0. However, we need to remem- ber that for the Gamma model, the prediction represents the mean of the total duration, while for the lognormal and mixed model it represents the median. Hence, we also calculate the mean of the total duration for the lognormal and mixed model via µ(D(RRUP= 0)) = exp µlnDS+ 0.5σP(RRUP= 0)2 (32) We calculate the mean for each posterior sample, which gives a posterior distribution for this parameter. In Figure 3, the posterior distribution of the mean prediction for RRUP= 0 are quite similar between the models. Figure 4 shows the estimated standard deviations/shape parameters. The mixed model shows a smaller value for σSthan the other models. The reason is that the event term δB in the mixed model is not just a source term, but models a different distance scaling for each event, and thus σSis more a measure of the variability of the distance scaling. This can be seen in Figure 5, which shows predictions for each of the 30 events (including the event-specific DSorδB). The Gamma and logormal (partially decoupled) models show a similar range of predictions, which stays constant with increasing distance. By contrast, the mixed model predictions show a different distance scaling for each event, and increasing range with distance. Within-event residuals are shown in Figure 6. The residuals for the Gamma and lognormal models are calculated via resid = lnDobs−ln [DS+fpath(⃗ r;RRUP)] (33) where DSis the estimated event-specific source duration. For the mixed model the residual is the term δWin Equation (28). Hence, the residuals are calculated in log-scale with respect to the logarithm of the mean prediction for the Gamma model, and with respect to the logarithm of the median prediction of the other two models. There are no large differences between the 10 01020 0.40 0.45 0.50 0.55 0.60 0.65r1densityGamma Lognormal Mixed 020406080 0.24 0.26 0.28 0.30r2density 0.000.030.060.09 50 60 70 80 90 Rbdensity 0.00.20.40.6 20 25 30 35 average source durationdensityGamma Lognormal Mixed Mean MedianFigure 3: Posterior density of distance slopes r1,r2, distance break point Rb, and average source duration, for the three different models. For the averge source duration, solid lines are the mean of the total duration at RRUP= 0, while dashed lines are the median. 0102030 0.0 0.1 0.2 0.3 0.4σSdensity 0.000.250.500.751.00 0 100 200 300 RRUP(km)σPLognormal Mixed 0510152025 0 100 200 300 RRUP(km)αPMean 5% fractile 95% fractile Figure 4: Posterior distribution of the standard Deviation σSof the source duration distribu- tion, as well as shape parameter αP(for the Gamma model) and scale parameter σP(for the lognormal and mixed model), dependent on distance. For the shape/scale parameters, the mean and 5%/95% fractiles are shown, calculated from the posterior distributions of the parameters associated with the shape/scale parameters. 050100150 0 100 200 300 RRUP(km)D5−75(s)Gamma 050100150 0 100 200 300 RRUP(km)D5−75(s)Lognormal 050100150 0 100 200 300 RRUP(km)D5−75(s)Mixed Figure 5: Individual mean (Gamma) and median (lognormal and mixed) prediction of the 30 different simulations for the three different models. Predictions are calculated with the mean of the posterior samples of the parameters. 11 −2−1012 0 100 200 RRUP(km)residualGamma −2−1012 0 100 200 RRUP(km)residualLognormal −2−1012 0 100 200 RRUP(km)residualMixedFigure 6: Residuals of the different models. The black line is a loess fit (Cleveland, 1979) to the residuals. 050100150200 0 100 200 300 RRUP(km)D5−75(s)Mean Median 5% fractile 95% fractile 31030100 10 30 100 300 RRUP(km)D5−75(s) Gamma Lognormal Mixed Figure 7: Predictive density and data for the different models. different models in terms of the residuals. The residuals of all models are well well centered, with larger variability at shorter distances. Such variability is also seen in peak and response spectral values in Frankel et al. (2018), and due to different distances to subevents. Figure 7 shows the distance scaling of the predictive distribution, together with the data. The predictive distribution comprises both the epistemic uncertainty and the aleatory vari- ability of the model. For each distance, we calculate the mean or median prediction for each posterior sample. Then, we randomly sample 100 source durations or event terms, and sample again 100 overall durations for each source duration/event term. From the resulting samples, we calculate the mean, median, and 5% and 95% fractiles. The mean and median predictions of all models are quite similar, but there are some differences in the fractiles, in particular at very small and very large distances. Figure 8 shows graphical posterior predictive checks (Gabry et al., 2019). For these plots, we draw samples from the posterior predictive distribution (cf. Equation (31)), and compare to the data. For each sample of the posterior distribution, we randomly draw a new data set (i.e. we calculate the mean/median and shape/standard deviation for each data point and sample from the corresponding distribution). Thus, we have 800 synthetic data sets of 15,810 simulated durations. In Figure 8, we compare the median, 5% and 95% fractiles across those 800 simulated data sets with the median and fractiles of the observed data. Overall, for a good model we would expect that the fractiles of the simulated data sets are similar to the ones data. In general, the differences between the models are small, but the Gamma model seems to be able to capture the 5% and 95% fractiles of the data slightly better. The plots shown in Figure 8 show median and fractiles of the whole data set(s). The models can perform differently for different ranges of the data; we provide plots similar to Figure 8 for different distance bins in the electronic supplement1. To quantify model comparison, we can look at the expected log pointwise predictive density 1available at https://github.com/nikuehn/Duration_M9 12 0.00.30.60.91.2 57 58 59 60 61 62 MediandensityGamma Lognormal Mixed 0.00.51.01.5 12 13 14 15 5% fractiledensity 0.00.10.20.30.40.5 120 125 130 95% fractiledensityFigure 8: Posterior predictive checks. Shown are densities of the median and 5%/95% fractiles of 800 simulated data sets. The black vertical line displays the median or fractile of the M9 data. Table 1: Differences in expected log pointwise predictive density ∆ delpdPSIS−LOO, based on PSIS-LOO, standard error of the difference, and difference in log-likelihood ∆ LLon test data when leaving out one event at a time. Model ∆delpdPSIS−LOO SE(∆delpdPSIS−LOO) ∆ LL SE(∆ LL) Gamma 0.0 0.0 0.0 0.0 Mixed -23.7 24.8 -52.0 9.9 Lognormal -242.9 20.3 -53.8 9.5 (elpd, Vehtari et al., 2017) for a new data set. Table 1 lists the differences in elpd for the different models, relative to the best model in terms of elpd. The elpd values are an approxi- mation to leave-one-out cross-validation. Thus, they are conditioned on the estimated source duration values/event terms. For details, see Vehtari et al. (2017), Vehtari et al. (2021), and the documentation of the loopackage (Vehtari et al., 2022)2. Table 1 shows that the Gamma model has a better predictive performance than the lognormal model, which is in line with the posterior predictive checks shown in Figure 8. The mixed model has similar elpdcompared to the Gamma model (within the standard error of the elpddifference), and performs much better than the lognormal model. Here one needs to remember that elpdis a measure of the predictive distribution conditional of the source durations or event terms, which in the case of the mixed model contain differences in the distance scaling for the events (Figure 5). The similar elpd values between the Gamma and mixed model, and the much larger value of the mixed model compared to the lognormal model indicates that modeling differences in distance scaling may be advantageous. To assess model performance on new events (i.e. with unknown source durations/event terms), we also perform a leave-one event-out test. Here, we estimate the models by leaving out all records from one event as test values, and then calculate the log-likelihood of the test data. The differences i the log-likelihood to the Gamma model are shown in Table 1, for five runs (i.e. five test events), calculated as ∆ LL=LLmod−LLGamma . For details on the calculation of the log-likelihood values of the different models, see the electronic supplement3. The Gamma model attains the highest log-likelihood value, but now the lognormal and mixed model are comparable, which is consistent with Figure 7. 2https://mc-stan.org/loo/index.html 3available at https://github.com/nikuehn/Duration_M9 13 5.1 NGA Sub Data The models estimated in this paper are estimated for 30 simulations with magnitude M= 9. Assuming a constant (magnitude independent) stress parameter ∆ σ, we can still make predictions for other magnitude values, via Equation (18). In this section, we calculate residuals to data from the NGA-Subduction (NGA-Sub) project (Bozorgnia et al., 2022; Contreras et al., 2022; Ahdi et al., 2022; Mazzoni et al., 2022, 2021). We use the interface data from the KBCG NGA-Sub model (Kuehn et al., 2020, 2022), and calculate residuals for the Gamma model (residuals to the other models are provided in the electronic supplement4). Predictions can be calculated as follows: the median stress parameter can be calculated from the estimated median source duration via µln ∆σ=−3 (µlnDS+ ln 4 .9 + 6 ln 10 + ln β) + (1 .5∗9 + 16 .05) ln 10 (34) which is just rearranging of Equation (18). For each record, we can then calculate the mean prediction and residuals as Dsource =µln ∆σ 101.5M+16.05−1/3 / 4.9 106β expσ2 2(35) µ=Dsource +fpath(⃗ r;R) (36) resid = lnY−lnµ (37) The residuals are shown in Figure 9 with respect to distance, magnitude and VS30. The residuals are clearly biased (the mean of the residuals is -0.887), so the model based on the M9 simulations is overpredicting the data. The distance scaling of the data at larger distances is similar to the model, but is quite different for short distances. Part of this could be due to the geometry of the subuction zone (the M9 simulations are done for Cascadia/Pacific Northwest, while the NGA-Sub data is mostly comprised of Japanese records). Another possibility is that the distance breakpoint RBis magnitude dependent, which would lead to a different distance scaling at short distances. Otherwise, the residuals show similar behavior as for the M9 simulations, with larger variability at short distances. The residuals do not show a great trend with magnitude, which might indicate that a constant stress drop is not an unreasonable model. The increase of the loess fit at large magnitudes is mainly due to the residuals of the Tohoku ( M= 9.12) event, which suggests a larger than average source duration for Tohoku (compared to the NGA Sub data). Durations have been found to be negatively correlated with PGA (Bradley, 2011) based on shallow active crustal ground motions. The event term of Tohoku is negative in Kuehn et al. (2020), which would imply a positive event term for duration, so the observation of Figure 9 is in line with expectation. The residuals show a positive trend at low values of VS30, indicating longer durations at softer sites. This trend is in agreement with previous studies (Bahrampouri et al., 2021; Kempton and Stewart, 2006; Afshari and Stewart, 2016). 6 Discussion We have presented and compared different parameterizations of additive duration models, based on the M9 simulations. The models differ in the placement of the random effect of the event term (affecting both source and path durations, as done previously in the literature, or only affecting the source duration), as well as in the observation likelihood (lognormal and Gamma). 4available at https://github.com/nikuehn/Duration_M9 14 −3−2−101 0 200 400 600 800 RRUP(km)residualGamma −3−2−101 5 6 7 8 9 MresidualGamma −3−2−101 100 300 1000 VS30(ms)residualGammaFigure 9: Residuals to Interface NGA-Sub data from Kuehn et al. (2020) with the Gamma model. The solid black lines show a loess fit (Cleveland, 1979) to the residuals, the dashed black lines are the mean of the residuals. We think that the partially decoupled model is more physical (or its interpretation is more in line with the physical intention). The reason is that the partially decoupled model does not introduce a different path scaling for each event, and thus the variability of the estimated source-specific effects represents source duration variability, and is not affected by variability in distance scaling. We find that the partially decoupled Gamma model performs best of the three models, both in terms of graphical posterior predictive checks (Figure 8) as well as when assessed by PSIS-LOO and leave-one event-out training training/test splits (Table 1). However, the differences between the model predictions are not large (Figure 7), which raises the question whether it matters which model is used. The mean and median predictions of all models are very similar, and differences in the fractiles of the predictive distribution only become pronounced a larger distances. Hence, while the partially decoupled Gamma model is preferable in terms of its physical interpretation and a slightly better predictive performance, this is only of concern of the full distribution of the model prediction is needed. Related to this point is that the predictive distribution is harder to compute for the partially decoupled models; in the mixed model, the total variance can be calculated as σ2 T=σ2 S+σ2 P, which is not valid for the decoupled models, and one has to resort to numerical or Monte Carlo integration. We also find that the Gamma model has better predictive performance compared to the lognormal model. Published duration models (both multiplicative and additive) generally model duration with a lognormal distribution. As stated before, we do not think that there is an a- priori reason to prefer one over the other. Both lognormally and Gamma distributed variables are positive, continuous, and right skewed, which makes these distributions reasonable choices for duration data. The two models are both similarly easy (or complex) to fit, with similar results, so we believe more research is necessary to assess which model is better. In particular, the Gamma model should be applied to a real data set of observed duration models, to check whether the conclusions drawn from the M9 simulations also apply there. There are advantages to using a lognormal model over the Gamma model. A lognormal model allows to easily model correlations of (log) duration residuals with residuals from PSA GMMs (Bradley, 2011). Such correlations can help to constrain duration model regressions (Stafford, 2008). In addition, the development of spatial correlation models for within-event residuals are easier with normally distributed residuals, though such models are not common (one exmple is Huang et al. (2020)). By contrast to the lognormal and mixed model, which model the median of the duration distribution, the Gamma model is parameterized by a function for the mean of the distribution. This requires some care when comparing models, but it is easy to convert from median to mean for the lognormal model. This raises a more general question: when a point estimate of ground- motion duration is used, which one should be used? Possible choices are the median, the mean, or the mode, which are all different fro the lognormal and Gamma distribution. Bora et al. 15 (2014) used the median prediction of their duration model to calculate response spectra via random vibration theory (RVT) in conjunction with an empirical Fourier amplitude spectrum (FAS) model, but switched to the mean in Bora et al. (2015) and Bora et al. (2019). Lavrentiadis and Abrahamson (2021) used the median duration ( D5−85of (Abrahamson and Silva, 1996) in their RVT calculations. Duration model predictions may be used in stochastic simulations (Boore, 1983, 2003), as for example in Sokolov and Chernov (2001) or (Bommer et al., 2017). Often it is not clear which duration (mean or median is used), but given that most models use a lognormal distribution, it is implicitly the median. Green et al. (2020) sample from the duration distribution to calculate liquefaction hazard. It is probably application dependent whether it makes sense to neglect duration variability (it might not matter, it might be computationally too demanding), but if the variability is taken into account, the different models can potentially lead to different results. As shown in Table 1, the predictive performance of the mixed model is much better compared to the lognormal model when conditioning on estimated event terms, but similar when testing against a new event. Since the event terms of the mixed model also account for event-specific distance scaling, this suggests the simplistic path model (with a constant distance scaling) can be improved. The 30 events are all simulated at the same stations, so the path effect is less of regional effect (like different scaling due to different paths because events are located in different regions), but is due to local differences in paths, which can arise due to different location, orientation and depth of the subevents used for simulation. Adding an event-specific distance scaling term to the partially decoupled models (modeled as a random effect to coefficient r2) leads to a large improvement with values of ∆ delpdPSIS−LOO= +265 .8 for the Gamma model and ∆delpdPSIS−LOO= +75 .1 for the lognormal model, compared to the Gamma model without event-specific distance scaling (the best model according to Table 1). When estimating a model on observed data such as the NGA-Sub data, we would expect differences in distance scaling between events/regions, similar to (partially) nonergodic effects in PSA models (Stafford, 2014). Several of the NGA-Sub models (e.g. Kuehn et al., 2022; Parker et al., 2021; Abrahamson and G¨ ulerce, 2022) include regional adjustment terms to account for differences in anelastic attenuation and linear site scaling. It is easy to include regional similar random effects structures in both the mixed and partially decoupled models. Moving to more fully nonergodic models (Anderson and Brune, 1999), incorporating a cell- specific distance scaling model (Kuehn et al., 2019; Dawood and Rodriguez-Marek, 2013) into both models should pose no conceptual problem. Nonergodic models based on spatially varying coefficient models (Bussas et al., 2017; Landwehr et al., 2016) can also be included, but should be done in the parameterization of Kuehn (2021), where the nonergodic latent effects are not integrated out. The models presented in this paper do not include site scaling, i.e. scaling whith VS30and/or Z1.0/2.5. In the literature of (mixed) additive duration models, site scaling has been included as an additive (Abrahamson and Silva, 1996; Kempton and Stewart, 2006) or multiplicative term (Bahrampouri et al., 2021; Afshari and Stewart, 2016). In general, the partially decou- pled models can accommodate both possibilities, and one should assess the implications when developing new models. We only briefly touched upon within-model epistemic uncertainty associated with the mod- els. The epistemic uncertainty is quantified via the posterior distribution of the parameters, which is shown in Figures 3 and 4. The uncertainty of the coefficients translates into uncertainty of model predictions, which generally should not be ignored. Plots of the epistemic uncertainty associated with the model predictions, similar to Figure 7, are shown in the electronic sup- plement (Figure S12)5. For the presented models, the epistemic uncertainty is rather small, 5available at https://github.com/nikuehn/Duration_M9 16 2468 0 100 200 300 RRUP(km)90% epistemic intervalGamma Lognormal Mixed 708090100110120 0 100 200 300 RRUP(km)90% predictive intervalFigure 10: 90% credible intervals for the epistemic model predictions and the full predictive distribution. which is due to the data distribution. There are 30 events, all of the same magnitude, and all records are simulated at the same stations, so the data coverage is quite extensive compared to an empirical model. The influence of the epistemic uncertainty on model predictions in the M9 duration models is small, but we would expect it to be larger for models based on observed data. Figure 10 shows the 90% credible intervals against distance, calculated as the 95% frac- tile minus the 5% fractile. Shown are the intervals for the epistemic uncertainty (Figure S12) and the full predictive distribution (Figure 7). At longer distances, the mixed and lognormal models show a larger difference in the predictive distribution compared to the Gamma model. The epistemic uncertainty of the mixed model increases stronger with distance than for the other models. This is due to the slightly wider posterior distribution of r2for the mixed model seen in Figure 3. As discussed before, the mixed model formulation leads to event-specific distance scaling, so we can interpret the r2coefficient as a mean distance scaling coefficient. It is informed by the distance scaling of the 30 events, and hence is based on less data, which leads to larger uncertainty. This is similar to the larger uncertainty of the global model in a partially nonergodic GMM compared to an ergodic global model (e.g. Kuehn et al., 2022). When applied to observed data, the M9-models are strongly biased (Figure 9), though the scaling with magnitude as well as the long-distance scaling seem to be in agreement with the data scaling. Hence, we do not advise to use the models directly, but we believe that it can be used to provide informative prior distributions for some parameters of an empirical subduction duration model. This paper is concerned with the estimation of additive duration models. We have illustrated the different models with the D5−75-significant duration. Some experiments indicate that similar conclusions can be drawn for D5−95, but we have not made a comprehensive analysis like for D5−75. We do not expect major differences for other significant duration measures apart from the values of the parameters. Absolute and bracketed durations (i.e. duration measures based on exceedance of a certain threshold acceleration) decrease with distance (e.g. Bommer et al., 2009; Stafford, 2008), so models for them require more adjustment. Code Availability The Stan and Rcodes to run the analyses are available at https://github.com/nikuehn/ Duration_M9 . 17 Acknowledgements Regressions are carried out using R-package cmdstanr (Gabry and ˇCeˇ snovar, 2021) and assessed with package posterior (B¨ urkner et al., 2022). ELPD values are calculated with package loo (Vehtari et al., 2022). Plots are made with ggplot2 (Wickham, 2016) and bayesplot (Gabry and Mahr, 2022). Loess fits are calculated with msir (Scrucca, 2011), and the tidyverse (Wickham et al., 2019) is used for data wrangling. We would like to thank Yousef Bozorgnia, Ken Campbell, Brian Chiou, Art Frankel, and the participants of the NGA-Subduction Project for helpful comments and discussions. We would also like to thank Art Frankel and Erin Wirth for providing the M9 durations. We would like to thank the participants of the NGA-subduction program for numerous discussions that led to improvements of the methodology presented in this paper. The NGA- Subduction project was supported by FM Global, the U.S. Geological Survey, the California Department of Transportation, and the Pacific Gas and Electric Company. The supports are gratefully acknowledged. The opinions, findings, conclusions or recommendations expressed in this publication are those of the authors and do not necessarily reflect the views of the study sponsors.