[orcid=0000-0003-3685-8124]

[1]

[1]

1]organization=Centro de Investigación en Matemáticas, addressline=Jalisco
S/N, city=Guanajuato, postcode=36023, country=Mexico

[orcid=0000-0002-5795-4345] [2] []

[3]

2]organization=Universidad Nacional Autónoma de México, city=Mexico City,
postcode=04510, country=Mexico

[cor1]Isaías Bañales

# Modelling Evolutionary Power Spectral Density Functions of Strong
Earthquakes Via Copulas

Isaías Bañales<EMAIL_ADDRESS>[  J. Andrés Christen<EMAIL_ADDRESS>Josué Tago<EMAIL_ADDRESS>[

###### Abstract

This paper proposes a new approach for analyzing seismic accelerograms using
the evolutionary Power Spectral Density function (ePSDF). The accelerogram of
an earthquake can be accurately modeled and simulated from its spectrogram,
based on the oscillatory stochastic processes theory. To adequately
characterize a spectrogram that is consistent with the response spectra, a
parametric model with 16 parameters is proposed. This model describes the
square of the amplitude spectrum, an envelope of the square of the
accelerogram, and a copula that constructs a time-frequency model from the
time and frequency marginals. The use of copulas to model a bivariate
probability distribution is a common practice in statistics, particularly when
the marginal distributions are known. The periodogram can be viewed as an
unnormalized probability density function, where the total energy serves as
the normalization constant, since the total energy of a seismic motion is
always finite. Additionally, a reduced model consisting of only 10 parameters
is presented, which may be especially valuable when only shear wave effects
are relevant.

###### keywords:

Spectral Density Oscillatory process Copula Stochastic simulation Response
Spectrum

##  1 Introduction

Analyzing earthquakes in the Cocos plate is a significant challenge for
Mexican seismology, since strong ground motions may cause severe structural
damage in Mexico City and other regions of the country. Several authors, such
as Arroyo et al. (2010), Ordaz et al. (1994) or Castro and Munguia (1993),
have made important advances in analyzing the amplitude spectra for the
southern region of Mexico.

As discussed in Kramer (1996), the response spectrum is a crucial ground
motion parameter in earthquake engineering. To assess seismic risk, it is
necessary to calculate response spectra directly from available earthquake
information.

One of the most influential papers for simulating accelerograms compatible
with response spectrum is Boore (2003); his algorithm, discussed in his
section Obtaining ground motions, works on Fourier domain after applying a
filter based on the envelope over time of the seismic acceleration. In that
paper a random vibration theory is proposed based on the main formula of
Cartwright and Longuet-Higgins (1956) which requires a stationary process, an
assumption that, evidently, earthquake accelerograms do not fulfill.

Additionally, certain decisions made in Boore (2003), such as the filtering
and windowing of the white noise, appear to lack clear justification. These
factors, combined with the assumption of stationarity, were the primary
motivations behind this study. We present a theoretical approach based on
oscillatory process (Priestley, 1965); the new approach allows us to work with
a theoretically well-based envelope of the earthquake over the time, in
contrast to the modulation in Boore (2003) based on "I find a good comparison
between response spectra computed using the box and exponential windows".

To avoid assuming a second-order stationary process, the oscillatory
stochastic processes theory is used as proposed by Priestley (1965). This
theory facilitates the utilization of spectrograms as effective estimators for
the Evolutionary Power Spectral Density Function (ePSDF) of an underlying
stochastic process responsible for generating the accelerogram. Moreover, the
spectrogram is used in this work to analyze the energy released in time and
frequency separately to finally connect them using the probabilistic theory of
copulas.

In the context of signal analysis, copulas have been used to model simulated
spectrograms of chirps, as presented in Boashash (2015) and Davy and Doucet
(2003).

A few novel papers have introduced the use of oscillatory stochastic processes
in seismology. In Vlachos et al. (2016) a mixture of Kanai-Tajimi filters with
time-dependent parameters is proposed, while in Broccardo and Dabaghi (2017) a
filter with frequency dependence is used.

In Vlachos et al. (2016) the idea is to sum multiple Kanai-Tajimi filters with
time-variant parameters. They present the next formula

$S(\mathfrak{f},t)=|\text{HP}(\mathfrak{f})|\sum_{k=1}^{2}S_{0}^{k}(t)\frac{1+\Big{(}2\xi^{k}(t)\mathfrak{f}/\mathfrak{f}^{k}_{g}(t)\Big{)}^{2}}{\Big{(}1-\Big{(}\mathfrak{f}/\mathfrak{f}^{k}_{g}(t)\Big{)}^{2}\Big{)}^{2}+\Big{(}2\xi^{k}(t)\mathfrak{f}/\mathfrak{f}^{k}_{g}(t)\Big{)}^{2}}$
(1)

where $\\{\xi^{k}(t),\mathfrak{f}^{k}_{g}(t),S_{0}^{k}(t)\\}$ are a time
varying modal dominant frequency, modal apparent damping ratio, modal
participation factor with respect to the $k$ mode, and HP is a high pass
filter.

The necessity of imposing a high pass filter in (1) arises from the fact that
Kanai-Tajimi filters do not have the value 0 at frequency 0. Also, the
modelling and estimation of
$\\{\xi^{k}(t),\mathfrak{f}^{k}_{g}(t),S_{0}^{k}(t)\\}$ produce a non
interpretable model that does not allow to incorporate previous models for the
amplitude frequency as the $\omega^{2}$ model for amplitude spectrum (Brune,
1970) or models for the accelerogram envelope as Stafford et al. (2009).

In the context of attenuation laws using ePSDF only the paper Vlachos et al.
(2018) is available. Its model proposes a 180 parameter regression. Also its
parameters are not interpretable, and they do not quantify the uncertainty of
their estimations.

The framework presented in this paper aims to simplify the attenuation law for
response spectra by first modeling the amplitude spectrum, the envelope of the
accelerogram over time, a copula, and the total energy of the accelerogram.
This approach results in an interpretable model that is informed by physical
considerations.

The seismic data were provided by the accelerographic network provided by the
Engineering Institute of UNAM (RAII UNAM, by its acronym in Spanish), product
of the instrumentation and processing work of the Seismic Instrumentation
Unit. The data is distributed through the Accelerographic Database System on
the web: https://aplicaciones.iingen.unam.mx/AcelerogramasRSM/.

The data used in section 4 is a subset of the records curated in García et al.
(2009) and used in Arroyo et al. (2010). The ground motions are interplate
thrust fault events with $M_{w}\geq 5$ and epicentral distance, $R$, between
100 and 500 km. For those records that did not contain $M_{w}$, they were
collected from the database of Global Centroid-Moment-Tensor (CMT) Project
(more information about the CMT project can be consulted in Ekström et al.
(2012)).

The outline of the paper is as follows. In section 2 the essential theory for
understanding the model is presented, which includes elementary results in
stochastic oscillatory processes and copulas and why they are useful to model
accelerograms. Section 3 presents a parametric way to describe an
accelerogram, as well as the physically informed constraints. Examples of the
model, a discussion of the errors, and a reduced version that works only with
the $S-$wave are presented in section 4. Section 5 talks briefly about further
work to incorporate attenuation laws into the framework presented in this
work. Finally section 6 summarizes the results obtained in this work.

##  2 Theoretical background

In this section, a brief review of the essential concepts for understanding
the model in section 3 is presented. The concept of oscillatory stochastic
processes is introduced, as well as a statistical method for estimating the
corresponding ePSDF. The ePSDF is a generalization of the classical spectral
density using Fourier analysis for non-stationary oscillatory processes.

###  2.1 Basic Theory on Oscillatory Stochastic Processes

In Priestley (1965) the concept of oscillatory processes is introduced for the
first time, which extends the family of stochastic processes admitting a
spectral representation.

If there is a complex stochastic process $\\{X(t):t\in\mathbb{R}\\}$, which is
trend free, i.e. $\mathbb{E}[X(t)]=0$ for all $t$, and if its covariance
function, defined by

$R_{s,t}=\mathbb{E}[X(s)X^{*}(t)],$

may be written as

$R_{s,t}=\int_{\mathbb{R}}\phi_{s}(\omega)\phi_{t}^{*}(\omega)d\mu(\omega),$
(2)

where $*$ denotes the complex conjugate, for functions
$\\{\phi_{t}(\omega):t\in\mathbb{R}\\}$ and a measure $\mu(\omega)$ defined on
the real line. Then, as it is presented in Priestley (1965), the process
admits the representation

$X(t)=\int_{\mathbb{R}}\phi_{t}(\omega)dZ(\omega),$ (3)

where $Z(\omega)$ is a self orthogonal process, i.e.

$Cov(Z(\omega),Z(\lambda))=0\quad\text{for all $\omega\neq\lambda$}$

and

$\mathbb{E}[(dZ(\omega))^{2}]=d\mu(\omega).$

In addition, in order to have $\text{Var}(X_{t})<\infty$ for all $t$, it is
required that $\phi_{t}$ is quadratically integrable. Over the paper it will
be assumed that $\text{Var}(X_{t})<\infty$.

If $\phi_{t}(\omega)=e^{i\omega t}$, the classical spectral representation of
an homogeneous stochastic process is retrieved. Further details about spectral
representations for stationary processes can be consulted in Wei (2006).

In Priestley (1965) the next definition is presented.

###### Definition 1

For a $t$ fixed, $\phi_{t}(\omega)$ will be said to be an oscillatory function
if for some (necessarily unique) $\theta(\omega)$, it may be written in the
form

$\phi_{t}(\omega)=A_{t}(\omega)e^{i\theta(\omega)t},$

where $A_{t}(\omega)$ is of the form

$A_{t}(\omega)=\int_{\mathbb{R}}e^{it\theta}dH_{\omega}(\theta),$

with $H_{\omega}(\theta)$ a measure that $\mid dH_{\omega}(\theta)\mid$ has an
absolute maximum at $\theta=0$.

If $\theta(\omega)$ is a single-valued function of $\omega$ then after a
variable change from $\omega$ to $\theta(\omega)$ in (3) and redefining
$A_{t}(w)$ and the measure $\mu(w)$, the process can be rewritten as

$X(t)=\int_{\mathbb{R}}A_{t}(\omega)e^{i\omega t}dZ(\omega).$

It is important to remark that $A_{t}(\omega)$ is a deterministic time-
frequency dependent function (as opposed to classic Fourier), which provides a
more flexible framework to model non-stationary processes.

###### Definition 2

If there exists a family of oscillatory functions $\\{\phi_{t}(\omega)\\}$, in
terms of which the process $\\{X(t)\\}$ has a representation of the form (3),
$\\{X_{t}\\}$ will be termed an oscillatory process.

In Priestley (1965) it is explained that if $0\leq t\leq T$ for $T$ in
$\mathbb{R}$, the representation (3) always exists. Considering this result,
the accelerograms can be assumed to be oscillatory processes. It is important
to note that this is only a sufficient condition, the theory of oscillatory
processes does not include only processes in a compact interval over the time.

In the case of an oscillatory process $\\{X_{t}\\}$ its variance can be
written as

$\textit{Var}(X(t))=R_{t,t}=\int_{\mathbb{R}}\mid
A_{t}(\omega)\mid^{2}d\mu(\omega),$

which can be seen as the analogous version of the representation in the
stationary case. The evolutionary power spectra $dF_{t}(\omega)$ is defined as

$dF_{t}(\omega)=\mid A_{t}(\omega)\mid^{2}d\mu(\omega),$

and if $\mu(\omega)$ is absolutely continuous respect to Lebesgue measure,
then

$dF_{t}(\omega)=f_{t}(\omega)d\omega,$

hence $f_{t}(\omega)$ will be referred to as the ePSDF.

###  2.2 Short Time Fourier Transform to Estimate the ePSDF

Since the representation (3) relates the processes $X(t)$ and
$\textit{Var}(X(t))$, the ePSDF can be used to simulate realizations of
$X(t)$.

A way to estimate the ePSDF is using the short time Fourier transform (STFT).
According to Allen (1977), the STFT is

$\text{STFT}(t,\omega)=\int w(t-\tau)X(\tau)e^{-i\omega\tau}d\tau,$ (4)

where $w$ is an appropriate window function. To guarantee that STFT preserve
the results of Theorem 7.2 in Priestley (1965), it is necessary that $w$
satisfies

$\int_{\mathbb{R}}w^{2}(t)dt=1,$ (5)

which guarantee that the estimator (7) is asymptotically almost unbiased and
also ensuring that the total energy of the signal does not change after
estimation. In this work $w(t)$ is of the form

$w(t)\propto e^{-\frac{t^{2}}{2(2k)^{2}}}\mathbbm{1}(\lvert t\rvert<k),$ (6)

for a $k$ in $\mathbb{R}^{+}$ previously defined, i.e. a Gaussian window and

$\mathbbm{1}(\lvert t\rvert<k)=\begin{cases}1&\text{if $\lvert t\rvert<k$}\\\
0&\text{if $\lvert t\rvert\geq k$}\\\ \end{cases}$

Based on Liang et al. (2007), the estimator of the ePSDF using a realization
$x_{o}(t)$ from $X(t)$ will be the square modulus of the STFT of $x(t)$, i.e.

$\overset{\bullet}{f}_{t}(\omega)=\Big{|}\int
w(t-\tau)x_{o}(\tau)e^{-i\omega\tau}d\tau\Big{|}^{2},$ (7)

also is defined

$\widetilde{f}_{t}(\mathfrak{f})=2\pi\overset{\bullet}{f}_{t}(\omega)$ (8)

where $\mathfrak{f}$ is expressed in Hertz instead of angular frequency.

To choose $k$ adequately in (6) it is necessary to understand the units of
$\widetilde{f}_{t}(\mathfrak{f})$, therefore in this section a dimensional
analysis is presented

The Parseval’s identity could be extended to oscillatory process, as is
presented in Liang et al. (2007) as

$E_{T}=\int_{0}^{\infty}x^{2}_{o}(t)dt=\int_{0}^{\infty}\int_{0}^{\infty}\widetilde{f}_{t}(\mathfrak{f})dtd\mathfrak{f},$
(9)

which ensures the conservation of total energy $E_{T}$ in our estimate
$\widetilde{f}_{t}(\mathfrak{f})$.

Since $X^{2}(t)$ has units in $\frac{cm^{2}}{s^{4}}$ (assuming $X(t)$ recorded
in $\frac{cm}{s^{2}}$), $E_{t}$ is in $\frac{cm^{2}}{s^{3}}$. Given this
consideration, we focus our interest in the units of
$\widetilde{f}_{t}(\omega)$ and its marginals. In (5) it is required that
$h^{2}(t)$ will be squared integrable and its integral will be 1, so $w(t)$
must be in $\frac{1}{\sqrt{s}}$, and by equation (4) it is easy to observe
that $\widetilde{f}_{t}(\mathfrak{f})$ is in $\frac{cm}{s^{3}}$, as expected
from (9). Having the units of $\tilde{f}_{t}(\mathfrak{f})$, it can be easily
asserted that the units of its marginals are $\frac{cm^{2}}{s^{4}}$ and
$\frac{cm^{2}}{s^{2}}$ respectively, where the marginals are defined as

$\displaystyle\tilde{g}_{1}(t)$
$\displaystyle=\int\tilde{f}_{t}(\mathfrak{f})d\mathfrak{f},$
$\displaystyle\tilde{g}_{2}(\mathfrak{f})$
$\displaystyle=\int\tilde{f}_{t}(\mathfrak{f})dt.$

where $\tilde{g}_{1}$ and $\tilde{g}_{2}$ denotes the marginals in frequency
and time respectively.

From equation (7) and considering an observed accelerogram, $x_{0}$, in the
time interval $(0,b)$, the value of $\tau$ where the integrated function is
non-zero is in $\Big{(}\max(\tau-k,0),\min(\tau+k,b)\Big{)}$, it can be
observed that when $k$ goes to infinite the expression in (6) turns into the
constant

$w(t)=\frac{1}{\sqrt{b}}$

which replaced in (8) results in

$\int_{0}^{b}\tilde{f}_{t}(\mathfrak{f})dt=\int_{0}^{b}\Big{|}\int_{\mathbb{R}}\frac{1}{\sqrt{b}}x_{o}(\tau)e^{-i2\pi\mathfrak{f}\tau}d\tau\Big{|}^{2}dt=A^{2}(\mathfrak{f}).$
(10)

For the marginal in time, when $k$ takes the value 0 it can be computed as

$\int_{\mathbb{R}^{+}}\tilde{f}_{t}(\mathfrak{f})d\mathfrak{f}=\int_{\mathbb{R}^{+}}\Big{|}\int_{\mathbb{R}}\delta_{t}(\tau)x_{o}(\tau)e^{-i2\pi\mathfrak{f}\tau}d\tau\Big{|}^{2}dt=x_{o}^{2}(t),$
(11)

where $\delta_{t}(\cdot)$ is a Dirac delta function in $t$.

In Allen (1977) it is discussed the requirement of a double application of the
Nyquist theorem to determine the number of times and frequencies needed to
represent the STFT. The ePSDF presents an uncertainty principle problem,
resulting in a trade-off between time and frequency resolution. Higher time
resolutions lead to poorer frequency resolutions and vice versa.

Taking into consideration the uncertainty principle, and equations (10) and
(11) a good option to choose the optimal time window, $k^{*}$, is

$k^{*}=\underset{k\text{ in }\mathbb{R}^{+}}{\mathrm{argmin}}\left\\{\|\lvert
x_{o}\rvert-\sqrt{\tilde{g}_{1}}\|_{2}^{2}+\|A-\sqrt{\tilde{g}_{2}}\|_{2}^{2}\right\\},$
(12)

where $\|\cdot\|_{2}^{2}$ is the squared of the $L_{2}$ norm.

According to Parseval’s theorem, each term in equation (12) represents the
energy loss of the signal generated by choice $k$ in both time and frequency
domains.

Following this approach and after experimenting with different earthquakes, a
time window of 4 seconds was chosen. To increase temporal precision,
overlapping intervals of $90\%$ were used. Vlachos et al. (2016) used time
windows of 3 seconds which produced similar results for the data used in this
work. An example of an ePSDF estimated from the STFT can be seen in Figure 1.

It should be noted that the marginal on time, $\tilde{g}_{1}(t)$, has squared
units with respect to the acceleration. On the other hand, the marginal on
frequency has squared units with respect to the amplitude spectrum ($A$). For
$k^{*}$ it can be seen in Figure 2 that $\sqrt{\tilde{g}_{1}(t)}$ and
$\sqrt{\tilde{g}_{2}(\mathfrak{f})}$ approximate very well the accelerogram
and the amplitude spectrum of the record respectively.

Figure 1:  Left: Accelerogram of the record. Right: relevant part of its
respective ePSDF as estimation of spectral density in Hz,
$\widetilde{f}_{t}(\mathfrak{f})$, where most of its density is concentrated.
For the record UNIO0301.221, of $M_{w}=7.5$ of 22 January 2003, recorded at La
Unión, state of Guerrero, Mexico, with an epicentral distance of $R=263.59km$
using the North-South chanel.

Figure 2: Left: $\sqrt{\tilde{g}_{1}(t)}$ over $x_{o}(t)$. Right:
$\sqrt{\tilde{g}_{2}(\mathfrak{f})}$ over $AS$. For the record UNIO0301.221

It is important to note that there are several approaches to estimating the
ePSDF besides the STFT, among them those based on wavelets or on the Hilbert-
Huang transform. In this work, a comparative evaluation of each possible
estimation strategy was unnecessary since Bruns (2004) proved that, under very
mild assumptions, the three approaches are the same. Other novel approaches
exist, Huang et al. (2022), Rosen et al. (2009) or Fraser and Boashash (1994),
but theorems demonstrating superior performance to STFT are still lacking. For
a review on this topic, the reader is referred to Abreu and Romero (2017).

###  2.3 Simulation of an Oscillatory Stochastic Process

The proposed model considers the accelerogram of a particular earthquake, at a
particular location, as a realization of a non-stationary stochastic process,
$X(t)$, the actually recorded observed accelerogram, $x_{o}(t)$, can be
considered as a realization of such stochastic process $X(t)$. Moreover, any
computer simulation $x(t)$ of the modeled accelerogram $X(t)$ is a different
realization of such stochastic process and may not be exactly equal to
$x_{o}(t)$. Note however, that it is not necessary to exactly recover
$x_{o}(t)$, since the model produces synthetic accelerograms, $x(t)$, that
have the same distributional properties properties as $X(t)$.

Conceptually, the approach is to view the recorded accelerograms as a
realization of an underlying stochastic process, and the aim is to describe
such an underlying stochastic process which, according to the previous
section, is fully described by its ePSDF $f_{t}(\mathfrak{f})$.

In Liang et al. (2007), an approach to the simulation of oscillatory processes
is presented. Using the representation given in (3), they observed that if
$X(t)$ is a real process, as it is in accelerograms, the conjugated process,
$X^{*}(t)$, must satisfy

$X^{*}(t)=X(t)=\int_{\mathbb{R}}A^{*}_{t}(\omega)e^{-i\omega
t}dZ^{*}(\omega)=\int_{\mathbb{R}}A^{*}_{t}(-\omega)e^{i\omega
t}dZ^{*}(-\omega).$

Therefore, $A$ and $dZ$ for all $\omega$ satisfy

$\displaystyle A_{t}^{*}(\omega)$ $\displaystyle=A_{t}^{*}(-\omega)$
$\displaystyle dZ(\omega)$ $\displaystyle=dZ^{*}(-\omega).$

Using this idea, Liang et al. (2007) prove that

$x_{o}(t)=\sqrt{2}\sum_{k=0}^{\infty}\Big{[}\overset{\bullet}{f}_{t}(\omega)\Delta\omega\Big{]}^{\frac{1}{2}}cos(\omega_{k}t+\Phi_{k}),$

where $w_{k}=k\Delta\omega$, $\Delta\omega\to 0$, and $\Phi_{k}$ are random
variables, independent and identically distributed as uniforms in $[0,2\pi]$.
In practice, it is used a finite sum as

$x(t)=\sqrt{2}\sum_{n=0}^{N}\Big{[}\overset{\bullet}{f}_{t}(\omega)\Delta\omega\Big{]}^{\frac{1}{2}}cos(\omega_{n}t+\Phi_{n}),$
(13)

where $\Delta\omega=\frac{\omega_{u}}{N}$, $\omega_{n}=n\Delta\omega$ and
$\omega_{u}$ is the Nyquist frequency. It is important to remember that the
evaluation must be performed at the times and frequencies required by the
model; in particular, the STFT has no continuity in time or frequency.

The Figure 3 shows an observed accelerogram is presented and three simulated
records generated from the equation (13) using $\widetilde{f}_{t}(\omega)$ of
an observed accelerogram.

Figure 3: Observed accelerogram for the same record presented in figure 1 and
three simulations using the estimated ePSDF presented in Figure 1 and,
simulating using equation (13).

From equation (13) it can be seen that the phase simulation is based solely on
the generation of uniform random variables, unlike Sato (2013) which proposes
the use of a fractional Brownian motion to simulate phases.

###  2.4 Copulas

Since $\tilde{f}_{t}(\mathfrak{f})$ is a non-negative function that preserves
the energy of $X(t)$, which is a finite number as presented in Liang et al.
(2007), it can be seen as a non-normalized probability density function
(p.d.f.). Accordingly, we define the p.d.f. version as

${f}_{t}(\mathfrak{f})=\frac{1}{Z}\widetilde{f}_{t}(\mathfrak{f}),$ (14)

where

$Z=\int_{0}^{\infty}\int_{0}^{\infty}\widetilde{f}_{t}(\mathfrak{f})dtd\mathfrak{f}.$

The advantage of working with $f_{t}(\mathfrak{f})$ as a p.d.f. is that
probability theory can be utilized, in particular the marginals in frequency
and time can be computed to recover the joint density using copulas. For a
review on copulas see Nelsen (2007). The definition of a copula and Sklar’s
theorem are also presented in the appendix A.

There are many copulas regularly used in applications (Nelsen, 2007). A widely
used copula is the Gaussian copula, which is associated with the Gaussian
distribution, and is defined as

$\displaystyle C_{\rho}(u,v)$
$\displaystyle=\mathcal{N}_{\rho}(\Phi^{-1}(u),\Phi^{-1}(v))$
$\displaystyle=\frac{1}{2\pi\sqrt{1-\rho^{2}}}\int_{-\infty}^{\Phi^{-1}(u)}\int_{-\infty}^{\Phi^{-1}(v)}\exp\Big{(}-\frac{1}{2}{\begin{bmatrix}s&t\end{bmatrix}\begin{bmatrix}1&\rho\\\
\rho&1\end{bmatrix}^{-1}{\begin{bmatrix}s\\\ t\end{bmatrix}}}\Big{)}dsdt,$

where $\mathcal{N}_{\rho}$ is the joint cumulative distribution function of a
bivariate Gaussian distribution with mean 0, variances 1, correlation $\rho$
and $\Phi^{-1}$ is the quantile function of a standard univariate Gaussian
distribution.

Section 3.3 uses a Gaussian copula to relate the time and frequency marginal
distributions. There are two important reasons for choosing a Gaussian copula.
First $\rho$ controls the correlation of the joint distribution and second
$C_{\rho}$ has a density with the following explicit formulation

$c_{\rho}(u,v)=\frac{1}{\sqrt{1-\rho^{2}}}\exp\Big{(}-\frac{1}{2}{\begin{bmatrix}\Phi^{-1}(u)&\Phi^{-1}(v)\end{bmatrix}\Big{(}\begin{bmatrix}1&\rho\\\
\rho&1\end{bmatrix}^{-1}-\begin{bmatrix}1&0\\\
0&1\end{bmatrix}\Big{)}{\begin{bmatrix}\Phi^{-1}(u)\\\
\Phi^{-1}(v)\end{bmatrix}}}\Big{)}.$

In Davy and Doucet (2003) there is a detailed discussion of the use of Copulas
for signal analysis. It is important to remark that Davy and Doucet (2003)
work with a complicated family of copulas with a complicated expression of the
density copula, which causes a difficult physical interpretation. To avoid
this, in this work it is proposed the Gaussian copula, which is less flexible,
but seems capable to model accelerograms correctly.

The Gaussian copula is useful for low-dimensional representation of
accelerogams because it has only one parameter, $\rho$. Section 3.3 discusses
how $\rho$ is restricted to create realistic models for accelerograms.
Additionally, the supplementary material presents accelerograms with different
conditions to emphasize that the Gaussian copula is flexible enough for our
purposes.

Moreover, other copulas were tested as the Clayton and Student’s t copulas
(Nelsen, 2007) as well as those presented in Chesneau (2023), among others. No
significant improvements where obtained with respect to the Gaussian copula,
to justify their more complex expressions.

##  3 A model to approximate the normalized ePSDF

This section presents an approach for modeling the marginal densities of a
normalized ePSDF on time and frequency. The strategy involves initially
modeling and approximating the marginals of the normalized ePSDF. Then, using
a Gaussian copula, an approximation of the full normalized ePSDF is obtained.
This, along with the calculated total energy, provides a distributional
description of the underlying stochastic process for the accelerogram. This
approximation will be parameterized with a total of 16 parameters, resulting
in a significant reduction in dimensionality compared to the full
accelerogram, which has $50-200$ samples per second.

###  3.1 Time domain marginal

To model the time marginal, it is required an envelope of the energy released
along time. In Stafford et al. (2009), a model based on a lognormal
distribution is presented to envelope a seismic wave. The lognormal
distribution $f_{[\mu,\sigma^{2},l]}$ is given by

$f_{[\mu,\sigma^{2},l]}(x)={\displaystyle{\frac{1}{(x-l)\sigma{\sqrt{2\pi}}}}\
\exp\left(-{\frac{\left(\ln\left((x-l)\right)-\mu\right)^{2}}{2\sigma^{2}}}\right)}\mathbbm{1}\\{x>l\\},$

the location parameter has been added to handle the problem of a non
standardized pre sample time in the records.

To take the $P$ and $S-$waves into account simultaneously, the time marginal
is model with a mixture of two lognormal distributions. To ensure that the $S$
wave releases more energy than the $P-$wave, phyisically informed restrictions
were imposed as is discussed in 3.4. The model for the time marginal is
formulated as

$\widehat{g}_{1}(t)=\pi_{1}f_{[\mu_{1},\sigma_{1}^{2},l_{1}]}(t)+(1-\pi_{1})f_{[\mu_{2},\sigma_{2}^{2},l_{1}+\gamma]}(t).$
(15)

The intention of the $\gamma$ parameter is to guarantee that the $S-$wave
arrives after the $P-$wave, namely $\gamma>0$. In Section 3.4, it is discussed
how prior information about wave velocities was used to further state
admissible values for $\gamma$. It is important to note that
$(\mu_{1},\sigma_{1}^{2},l_{1})$ is related to the $P-$wave and
$(\mu_{2},\sigma_{2}^{2},l_{1}+\gamma)$ to the $S-$wave.

A further extension of this work could be the addition of prior information
based on partition energy among elastic waves as in Sánchez-Sesma et al.
(2011).

###  3.2 Frequency domain

A seminal paper on amplitude spectrum modelling in a Bayesian framework in
Mexico is Ordaz et al. (1994), which uses the $\omega^{2}$ model. A brief
survey of amplitude spectrum models is presented in Boore (2003) and a model
for finite-source in Mexico can be consulted in Singh et al. (1989).

The amplitud frequency models presented in Boore (2003) could be well
approximated by a lognormal distribution. To model the $P-$wave frequency
marginal a single lognormal distribution is used, while for the $S$ case a
mixture of two lognormal distributions is proposed. The physical reason to do
this distinction is that $S-$waves could be overlapped producing a multimodal
$S-$wave marginal or because the $S-$wave could transport higher energy on a
specific frequency. The models for the frequency marginals are

$\displaystyle\widehat{g}_{21}(\mathfrak{f})$ $\displaystyle\propto
f_{[\mu_{\mathfrak{f}_{1}},\sigma_{\mathfrak{f}_{1}}^{2},0]}(\mathfrak{f})\mathbbm{1}_{\mathfrak{f}\leq\mathfrak{f}_{max}}$
$\displaystyle\widehat{g}_{22}(\mathfrak{f})$
$\displaystyle\propto\Big{(}m_{2}f_{[\mu_{\mathfrak{f}_{21}},\sigma_{\mathfrak{f}_{21}}^{2},0]}(\mathfrak{f})+(1-m_{2})f_{\mu_{[\mathfrak{f}_{22}},\sigma_{\mathfrak{f}_{22}}^{2},0]}(\mathfrak{f})\Big{)}\mathbbm{1}_{\mathfrak{f}\leq\mathfrak{f}_{max}},$

The $\widehat{g}_{21}$ is the model for the frequency marginal before the
$S-$wave arrival time, i.e. before $l_{1}+\gamma$, and $\widehat{g}_{22}$
after the arrival. Note that the location parameters are fixed at 0 unlike as
in the time marginal density.

As presented in Boore (2003), Arroyo et al. (2010) and Parajuli (2017), it was
chosen $\mathfrak{f}_{max}=15$Hz as higher frequencies are usually not of
interest.

###  3.3 Time-Frequency domain model using copulas

As explained earlier, with the approximations for the time and frequency
marginals, the joint density is approximated using a Gaussian copula.

The proposal is to use a Gaussian copula with negative correlation to recover
the joint distribution $f_{t}(\mathfrak{f})$, with the following model, namely

$\displaystyle\widehat{f}_{t}(\mathfrak{f})$
$\displaystyle:=\widehat{G}_{1}(l_{1}+\gamma)c_{\rho}(\widehat{G}_{11}(t),\widehat{G}_{21}(\mathfrak{f}))\widehat{g}_{11}(t)\widehat{g}_{21}(\mathfrak{f})$
(16)
$\displaystyle+(1-\widehat{G}_{1}(l_{1}+\gamma))c_{\rho}(\widehat{G}_{12}(t),\widehat{G}_{22}(\mathfrak{f}))\widehat{g}_{12}(t)\widehat{g}_{22}(\mathfrak{f}),$

The densities $g_{1i}$ for $i$ in $\\{1,2\\}$ are defined as

$\displaystyle\widehat{g}_{11}(t)\propto\widehat{g}_{1}(t)\mathbbm{1}(t<l_{1}+\gamma)$
$\displaystyle\widehat{g}_{12}(t)\propto\widehat{g}_{1}(t)\mathbbm{1}(t>l_{1}+\gamma),$

to be probability density functions. The capital letters are the distribution
functions of the corresponding density functions, for example

$\widehat{G}_{1}(x)=\int_{0}^{x}\widehat{g}_{1}(u)du.$

In Section 4, there are examples of the fitting of the proposed model to real
seismic motions.

It is important to note that in order to calculate the response spectra, it is
necessary to know the energy of the record, which is easily calculated from

$E_{T}=\int X^{2}(t)dt.$ (17)

###  3.4 Model Fitting

Equation (16) requires the estimation of following parameters

$\Theta=\\{\mu_{1},\sigma^{2}_{1},l_{1},\mu_{2},\sigma_{2}^{2},\gamma,\pi_{1},\mu_{\mathfrak{f}_{1}},\sigma^{2}_{\mathfrak{f}_{1}},\mu_{\mathfrak{f}_{21}},\sigma^{2}_{\mathfrak{f}_{21}},\mu_{\mathfrak{f}_{22}},\sigma^{2}_{\mathfrak{f}_{22}},\pi_{2},\rho\\}.$

Since the signals have different orders of magnitude before and after the time
$l_{1}+\gamma$, the fitting of the parameters were accomplished using additive
Gaussian errors with different variance in each block to define a likelihood
in which the variance of the error is also estimated.

Defining $\mathcal{T}$ and $\mathcal{F}$ as the values in which the time and
frequency can be evaluated in the STFT, the negative log-likelihood is given
by

$\displaystyle-\ell(\Theta)$
$\displaystyle=\sum_{t\in\mathcal{T}}\sum_{\mathfrak{f}\in\mathcal{F}}\frac{\Big{(}\widehat{f}_{t}(\mathfrak{f})-f_{t}(\mathfrak{f})\Big{)}^{2}}{2\sigma_{\epsilon_{1}}^{2}}\mathbbm{1}(t<l_{1}+\gamma)\mathbbm{1}(\mathfrak{f}<\mathfrak{f}_{max})$
(18)
$\displaystyle+\sum_{t\in\mathcal{T}}\sum_{\mathfrak{f}\in\mathcal{F}}\frac{\Big{(}\widehat{f}_{t}(\mathfrak{f})-f_{t}(\mathfrak{f})\Big{)}^{2}}{2\sigma_{\epsilon_{2}}^{2}}\mathbbm{1}(t\geq
l_{1}+\gamma)\mathbbm{1}(\mathfrak{f}<\mathfrak{f}_{max})$
$\displaystyle+\Big{(}\sum_{\mathfrak{f}\in\mathcal{F}}\mathbbm{1}(\mathfrak{f}<\mathfrak{f}_{max})\Big{)}\Big{(}\sum_{t\in\mathcal{T}}\log(\sigma_{\epsilon_{1}}^{2})\mathbbm{1}(t<l_{1}+\gamma)$
$\displaystyle+\sum_{t\in\mathcal{T}}\log(\sigma_{\epsilon_{2}}^{2})\mathbbm{1}(t\geq
l_{1}+\gamma)\Big{)}$

To compute a synthetic accelerogram from the model, it is necessary to know
the total energy, $E_{t}$, which can be computed directly from the
acceleration record, $X(t)$ (see equation (17)).

It can be observed that if the Gaussian error $\sigma_{\epsilon_{i}}$ in each
section is known, the equation (18) produces a weighted least squares fitting.

It is important to note that $\sigma_{\epsilon_{1}}$ and
$\sigma_{\epsilon_{2}}$ are nuisance parameters, so the model needs only 16
parameters to describe the entire history of an accelerogram, which are

$\\{\mu_{1},\sigma^{2}_{1},l_{1},\mu_{2},\sigma_{2}^{2},\gamma,\pi_{1},\mu_{\mathfrak{f}_{1}},\sigma^{2}_{\mathfrak{f}_{1}},\mu_{\mathfrak{f}_{21}},\sigma^{2}_{\mathfrak{f}_{21}},\mu_{\mathfrak{f}_{22}},\sigma^{2}_{\mathfrak{f}_{22}},\pi_{2},\rho,E_{T}\\}.$

To guarantee that $\pi_{i}$ in $i=1,2$ are a mixture parameters, they are
constrained over the interval $(0,1)$.

In order to have a negative correlation in the Gaussian copula $\rho$ is
constrained to be less than 0. The idea that supports the use of a negative
correlation between time and frequency is that the Earth is an anelastic
medium that acts as a low-pass filter, resulting in that lower frequencies
prevail more time than the higher frequencies. A detailed discussion regarding
energy release in seismic waves can be found in Shearer (2009).

The velocity model of southern Mexico was taken from Castro and Munguia (1993)
with $P-$wave velocity of $6.1\frac{km}{s}$ and $S-$wave velocity of
$3.5\frac{km}{s}$.

The arrival times can be roughly estimated using epicentral distance, $R$, and
the corresponding wave velocity. The difference between the estimated arrival
times of the $P$ and $S-$waves

$T_{b}=\frac{R}{3.5}-\frac{R}{6.1}$

is used to constrain $\gamma$ in $(0.5T_{b},1.5T_{b})$, letting the actual
arrival time difference (i.e. $\gamma$) to be within 50% above and below the
simple estimate $T_{b}$.

To avoid unidentifiability and to guarantee that the first wave corresponds to
the $P-$wave, it is forced in (15) that the density
$f_{\mu_{1},\sigma^{2}_{1},l_{1}}$ reaches its maximum before
$f_{\mu_{2},\sigma^{2}_{2},l_{1}+\gamma}$ does. Another constraint is that the
peak of $f_{\mu_{1},\sigma^{2}_{1},l_{1}}$ and
$f_{\mu_{2},\sigma^{2}_{2},l_{1}+\gamma}$ must be in the interval
$(0.5T_{b},1.5T_{b})$, to guarantee that the first energy arrival occurs
before the second.

Finally, to avoid too flat densities, the variances of the two lognormals are
bounded by $50^{2}$. This bound is arbitrary and can be improved if a model of
the duration of seismic waves is used. Since the $P-$wave register contains
less energy, the last condition
$0<\sigma_{\varepsilon_{1}}^{2}<\sigma_{\varepsilon_{2}}^{2}$ is imposed.

Then the optimization problem to obtain $\Theta$ is

$\displaystyle\underset{\Theta}{\mathrm{argmin}}-\ell(\Theta)$ subject to
$\displaystyle\Theta>0,$ $\displaystyle\pi_{i}\text{ in }(0,1),\text{ for $i$
in }\\{1,2\\},$ $\displaystyle\sigma_{\epsilon_{1}}<\sigma_{\epsilon_{2}},$
$\displaystyle\gamma\text{ in }(0.5T_{b},1.5T_{b}),$
$\displaystyle\max_{t}(f_{\mu_{1},\sigma^{2}_{1},l_{1}}(t))<\max_{t}(f_{\mu_{2},\sigma^{2}_{2},l_{1}+\gamma}(t)),$
$\displaystyle\underset{t}{\mathrm{argmax}}(f_{\mu_{2},\sigma^{2}_{2},l_{1}+\gamma}(t))-\underset{t}{\mathrm{argmax}}(f_{\mu_{1},\sigma^{2}_{1},l_{1}}(t))\text{
in }(0.5T_{b},1.5T_{b})\text{ and }$
$\displaystyle{\displaystyle[\exp(\sigma_{i}^{2})-1]\exp(2\mu_{i}+\sigma_{i}^{2})}<50^{2}\text{
for $i$ in }\\{1,2\\},$

where $\Theta>0$ means that each element of vector $\Theta$ is positive.

This is a difficult nonlinear optimization problem with complex constraints.
For the examples presented here the t-walk algorithm (Christen and Fox, 2010)
was used to carry out the required optimization . The t-walk is actually a
Markov Chain Monte Carlo (MCMC) algorithm, but it can easily be used for
optimization as well. It is derivative free, does not require tuning
parameters, domain restrictions are straightforward to include and has several
implementations, including one in Python.

##  4 Examples

In this section, the fitting of the model to real records will be shown. Only
two examples are presented in detail but several more can be found in the
supplementary material.

###  4.1 Full model examples

Figures 4 and 6 show the fit of the proposed model to a horizontal component
of the records ATYC8509.191 and UNIO0301.221.

The fitted model based on the copula approach presents satisfactory results
for marginal densities on time and frequency, which means that the energy of
our model is correctly released. Also the response spectra of the simulated
records approximates quite well to the real response spectra, as seen in
Figures 5 and 7.

Based on the different examples presented in this section and in the
supplementary material, the proposed model is flexible and useful for modeling
different seismic waves.

One of the major advantages of having a model to simulate accelerograms is
that it allows the analysis of how different wave paths change the response
spectra at a specific station, even for events with similar $M_{w}$ and $R$.
The response spectra of the records shown in figures 4 and 6 may be seen in
figures 5 and 7, respectively.

Figure 4: Top left: A simulated accelerogram using our fitted model (orange)
plotted over the observed accelerogram (blue). Top right: Time marginal
$g_{1}(t)$ and $\hat{g}_{1}(t)$. Bottom left: Frequency marginal of
$g_{2}(\mathfrak{f})$ and $\hat{g}_{2}(\mathfrak{f})$. Bottom right:
$f_{t}(\mathfrak{f})$. For the record ATYC8509.191, of $M_{w}=7.5$ of 19
September 1985, recorded at Atoyac, state of Guerrero, Mexico, with an
epicentral distance of $R=283.09km$ using the North-South chanel.

Figure 5: Same accelerogram as in Figure 4. Left: level curves of our fitted
full model overlapped to the ePSDF. Right: the response spectra of 20
simulations using our fitted model, and their median, and the response spectra
of the observed accelerogram.

Figure 6: Top left: A simulated accelerogram using our fitted model (orange)
plotted over the observed accelerogram (blue). Top right: Time marginal
$g_{1}(t)$ and $\hat{g}_{1}(t)$. Bottom left: Frequency marginal of
$g_{2}(\mathfrak{f})$ and $\hat{g}_{2}(\mathfrak{f})$. Bottom right:
$f_{t}(\mathfrak{f})$. For the record UNIO0301.221

Figure 7: Same accelerogram as in Figure 6. Left: level curves of our fitted
full model overlapped to the ePSDF. Right: the response spectra of 20
simulations using our fitted model, and their median, and the response spectra
of the observed accelerogram.

###  4.2 Reduced model

The response spectrum analysis usually takes only the $S-$wave and ignore the
$P-$wave Ordaz et al. (1994), Arroyo et al. (2010), Vlachos et al. (2016),
Tang et al. (2016), Sato (2013). To consider only the $S-$wave, the model can
be simplified in a straightforward way. To work only with the $S-$wave, ignore
the values before $l_{1}+\gamma$ and describe the accelerogram only with the
$S-$wave parameters. Namely

$\Theta^{*}=\\{\mu_{2},\sigma_{2}^{2},\pi_{1},\mu_{\mathfrak{f}_{21}},\sigma^{2}_{\mathfrak{f}_{21}},\mu_{\mathfrak{f}_{22}},\sigma^{2}_{\mathfrak{f}_{22}},\pi_{2},\rho,E_{T}\\},$

reducing the number of parameters to only 10. This is considerably less than
the 18 parameters in the $S-$wave modeling approach of Vlachos et al. (2016).

Figures 8 and 10 present examples of the reduced model with the same
accelerograms presented in Figures 4 and 6, respectively. It is important to
note that although the reduced model is less flexible, the response spectra,
in Figures 9 and 11, are comparable with those obtained from the full model,
in Figures 5 and 7 .

There is a trade off to take into account, the reduced model is faster and
easier to fit and allows for more parsimonious attenuation laws, but it
describes the ePSDF less accurately, which will slightly increase the error.
Depending on the importance of the $P-$available computing power, the reduced
model may be preferred over the full model.

Figure 8: Summary of the fit obtained for our reduced model for the 19
September 1985 earthquake, with $M_{w}=8.0$ and $R=283.09km$ recorded at
Atoyac, Guerrero. The observed data is presented in blue, the simulation and
the fit is presented in orange.

Figure 9: Left:level curves of the reduced model overlapped to the ePSDF.
Right: Response spectra of 20 distinct simulations using the reduced model, in
black its median and in blue the response spectra of the real accelerogram for
the record of Figure 8

Figure 10: Summary of the fit obtained for our reduced model for the 22
January 2003 earthquake, with $M_{w}=7.5$ and $R=263.59km$ recorded at La
Unión, Guerrero. The observed data is presented in blue, the simulation and
the fitting is presented in orange.

Figure 11: Left:level curves of the reduced model overlapped to the ePSDF.
Right: Response spectra of 20 distinct simulations using the reduced model, in
black its median and in blue the response spectra of the real accelerogram for
the record of Figure 10

###  4.3 Error analysis

To analyze the performance of the full model, the response spectra using the
synthetic accelerograms are compared to the observed ones. The procedure was
as follows. For each recorded accelerogram, $x_{o}$, 20 synthetic
accelerograms, $x_{i}$ for $i=1,2,...,20$, were simulated using
$E_{T}\hat{f}_{t}(\mathfrak{f})$ as their ePSDF. Subsequently, with the
synthetic accelerograms, the response spectra,
$\\{S_{a}(\omega_{0},0.05,x_{i}):i=1,2,...,20\\}$, were evaluated for the
natural periods of resonance $T$ in $\\{0.5,1,1.5,2,3,4\\}$. Finally, the
median of the response spectra at time $T$ for the 20 simulations,
$m(T,\\{x_{i}:i=1,...,,20\\})$, was compared with the value obtained from the
observed accelerogram.

Figure 12 displays the errors of the proposed model for each record in our
database. The error is simply defined as

$E(\cdot,x_{o})=S_{a}(\cdot,0.05,x_{o})-m(\cdot,\\{x_{i}:i=1,2,...,20\\}).$

The errors do not have trend and almost all of them are less than $20$ gals.
Moreover, the error band decreases with the period, reaching a band of $5$
gals for a period of $4$ s. These results are competitive with the errors in
the response spectra shown in Arroyo et al. (2010) and Vlachos et al. (2016).

It is important to note that the proposed model assumes a parametric
representation of a non-stationary process, which has fewer parameters than
similar ePSDF models, and that the inclusion of a priori physical knowledge in
the marginals allows for a flexible and competitive model with parameters that
are easy to interpret.

Reviewing the errors in the figure 12, the largest error corresponds to the
record UNIO8509.191 for the period of $0.5$ s with a value of $-64.90$cm/s2.
When the relative error is calculated, it turns out that the estimate varies
only 21% with respect to the response spectrum data. The error is large
compared to other earthquakes due to the high accelerations recorded for an
earthquake of $M_{w}=8$. In the supplementary material a good approximation of
the ePSDF can be seen, although the energy release was highly scattered.
Similar situations occur for the records PAPN8509.191 and SUCH8509.19 of the
19 September 1985 earthquake.

The records ATYC8509.211 and COYC8509.211 are from the earthquake of 21
september 1985 earthquake, which had a magnitude of $M_{w}=7.5$. In both
cases, our model underestimated the frequencies below 10 Hz, resulting in a
smaller predictions of the response spectra.

(a) Period of $0.5s$.

(b) Period of $1s$.

(c) Period of $1.5s$.

(d) Period of $2s$.

(e) Period of $3s$.

(f) Period of $4s$.

Figure 12: Difference of the response spectra of the observed accelerogram
versus the median of the response spectra of 20 simulation using the full
model.

##  5 Further work on attenuation laws

Each pair of $(\mu_{i},\sigma_{i}^{2})$ of a lognormal distribution is
information about the location and variability of the energy released, which
allows us to easily understand the physical meaning of our parameters in
$\Theta^{*}$, which is an important advantage over the model parameterization
of Vlachos et al. (2016) or the non-parametric model in Broccardo and Dabaghi
(2017).

To use our model in the context of attenuation laws, it will be necessary to
predict $E_{T}$ from $M_{w}$ and $R$. As a first approach it is proposed the
heuristic relation

$E_{T}=M_{w}e^{0.1R}.$ (19)

Figure 13(a) is a scatter plot of our proposed relation, showing that simple
regression models could predict $E_{T}$.

For the lag time between waves, Figure 13(b) is a scatter plot of the
epicenter distance versus lag time between copulas, i.e. the $\gamma$
parameter, which shows a useful relationship for predicting $\gamma$ from $R$
only.

With this evidences, it is encouraging to further extend our model for
attenuation laws in the context of seismic risk, since this novel approach for
the attenuation laws reach an important complexity reduction compared with
Vlachos et al. (2018) and it avoids the assumption of stationarity in contrast
to the papers based in the methodology of Boore (2003)

(a) $(E_{T},M_{w}e^{0.1R})$.

(b) $(\log(R),\log(\gamma))$.

Figure 13: The scatter plots were made using all the records.

##  6 Conclusions

The development of a fully physically inspired model that considers $P$ and
$S$ waves simultaneously using only 16 scalar parameters to analyze an
accelerogram was presented. Furthermore, the amplitude spectrum or seismic
envelope models could be used to generate attenuation models. The framework
presented also provides a flexible and easily interpretable model that
produces errors that are competitive with other recently developed models.

The main advantage of the presented framework is that, given the copula, the
time and frequency marginals are independent, allowing each marginal to be
modeled separately to finally obtain the joint distribution with a copula.

Unlike the models based on Kanai-Tajimi filters, which require an additional
high-pass filter to reduce the energy content in the lower frequencies, the
model presented in this paper allows the energy decay for the lower
frequencies to be controlled by the probability density function.

Even though the knowledge of the marginals in our model is designed for the
specific region of the Southern Mexico, the model is flexible enough to be
easily modified to work with other regions.

Implementing these ideas for ground motions with $R>500$, the addition of
surface waves could easily be done by adding more elements to the mixtures of
the marginals or using probability density functions with light tails for the
frequency margin to concentrate the energy in a smaller range of frequencies.

Although our model uses the $P$ and $S$ wave arrival times to constrain the
optimization problem, if they are known, they could be fixed for a better
inference. It is also important to note that the full model does not try to
identify each wave with high accuracy, our model only uses the information of
the waves to create an adequate envelope of the earthquake.

All the codes presented in this work can be found in the repository
https://github.com/isaiasmanuel/Attenuation_Laws and are available for
download.

##  Appendix A Basics of Copulas

In order to have a self-contained text, the definition of a copula and the
Sklar’s theorem are presented.

###### Definition 3

A two-dimensional copula is a function $C$ with the following properties:

  1. 1.

$C$ domain is $[0,1]\times[0,1]$

  2. 2.

$C(x,0)=0=C(0,y)$ for all $(x,y)$ in $[0,1]\times[0,1]$

  3. 3.

If $0\leq x_{1}\leq x_{2}\leq 1$ and $0\leq y_{1}\leq y_{2}\leq 1$ then

$C(x_{2},y_{2})-C(x_{2},y_{1})-C(x_{1},y_{2})+C(x_{1},y_{1})\geq 0$

  4. 4.

For every $(u,v)$ in $[0,1]\times[0,1]$ is met

$C(u,1)=u\quad\text{and}\quad C(1,v)=v.$

A copula is usually interpreted as a bivariate distribution function on the
unitary square with uniforms in $[0,1]$ as marginal distributions.

A key theorem in copulas theory is Sklar’s theorem, which relates every
bivariate distribution function to an associated copula; moreover, to define a
bivariate continuous distribution, it is only necessary to know the marginal
distributions and the associated copula.

###### Theorem 1

Sklar’s Theorem. Let $H$ be a joint distribution function with margins $F$ and
$G$. Then there exists a copula $C$ such that for all $x,y$ in
$\overline{\mathbb{R}}$,

$H(x,y)=C(F(x),G(y)).$ (20)

If $F$ and $G$ are continuous, then $C$ is unique; otherwise, $C$ is uniquely
determined on the cross product $\text{rank}(F)\times\text{rank}(G)$.
Conversely, if $C$ is a copula and $F$ and $G$ are distribution functions,
then the function $H$ defined by (20) is a joint distribution function with
marginal distributions $F$ and $G$.

The corollary 1 gives a direct way to recover the copula associated to a joint
distribution, it is very useful to define copulas with properties associated
to a specific distribution.

###### Corollary 1

Let $H,F,G$ and $C$ as in Theorem 1 and let $F^{(-1)}$ and $G^{(-1)}$ be the
quantile functions of $F$ and $G$, respectively. Then for any $(u,v)$ in
$[0,1]\times[0,1]$,

$C(u,v)=H(F^{-1}(u),G^{-1}(v)).$

##  Appendix B Response spectrum

A very important ground motion parameter is the response spectrum, which is
discussed in detail in Kramer (1996). A main idea during the development of
the model was to describe seismic motion that can reproduce the response
spectrum. In this section, the response spectrum is briefly introduced.

For an angular frequency $\omega_{0}$ and a damping ratio $\zeta$, the
displacement $u(t)$ is modeled by a forced damped harmonic oscillator, i.e. it
follows

$\frac{d^{2}u}{dt^{2}}+2\zeta\omega_{0}\frac{du}{dt}+\omega_{0}^{2}u=-x(t),$
(21)

with initial conditions

$\displaystyle u(0)=0$ $\displaystyle\frac{d}{dt}u(0)=0,$

where $x$ is an accelerogram.

The response acceleration spectrum of $x$ is
$\max_{t}\mid\frac{d^{2}u}{dt^{2}}\mid.$

As it is discussed in Kramer (1996) is common to take the pseudo acceleration
response spectrum (PARS), defined by

$S_{a}(\omega_{0},\zeta,x)=\omega_{0}^{2}\max_{t}\mid u\mid,$ (22)

which is an approximation of the response acceleration spectrum. In this work
the response spectra were computed with equation (22).

To solve the equation (21) the SciPy extension LSODA on Python was used, in
where Adams and BDF method with automatic stiffness detection and switching
are implemented. The damping parameter was set to $\xi=0.05$, as it is the
common value used in the specialized literature, as in Arroyo et al. (2010),
Boore (2003), Bora et al. (2016), Barone et al. (2015), Alimoradi and Beck
(2015), among others.

## References

  * Abreu and Romero (2017) Abreu, L.D., Romero, J.L., 2017\.  Mse estimates for multitaper spectral estimation and off-grid compressive sensing.  IEEE Transactions on Information Theory 63, 7770–7776.
  * Alimoradi and Beck (2015) Alimoradi, A., Beck, J.L., 2015\.  Machine-learning methods for earthquake ground motion analysis and simulation.  Journal of Engineering Mechanics 141, 04014147.
  * Allen (1977) Allen, J., 1977.  Short term spectral analysis, synthesis, and modification by discrete fourier transform.  IEEE Transactions on Acoustics, Speech, and Signal Processing 25, 235–238.
  * Arroyo et al. (2010) Arroyo, D., García, D., Ordaz, M., Mora, M.A., Singh, S.K., 2010.  Strong ground-motion relations for mexican interplate earthquakes.  Journal of Seismology 14, 769–785.
  * Barone et al. (2015) Barone, G., Lo Iacono, F., Navarra, G., Palmeri, A., 2015\.  A novel analytical model of power spectral density function coherent with earthquake response spectra, in: UNCECOMP 2015 - 1st ECCOMAS Thematic Conference on Uncertainty Quantification in Computational Sciences and Engineering, pp. 394–406.  doi:10.7712/120215.4280.805.
  * Boashash (2015) Boashash, B., 2015.  Time-frequency signal analysis and processing: a comprehensive reference.  Academic press.
  * Boore (2003) Boore, D.M., 2003.  Simulation of ground motion using the stochastic method.  Pure and applied geophysics 160, 635–676.
  * Bora et al. (2016) Bora, S.S., Scherbaum, F., Kuehn, N., Stafford, P., 2016\.  On the relationship between fourier and response spectra: Implications for the adjustment of empirical ground-motion prediction equations (gmpes).  Bulletin of the Seismological Society of America 106, 1235–1253.
  * Broccardo and Dabaghi (2017) Broccardo, M., Dabaghi, M., 2017\.  A spectral-based stochastic ground motion model with a non-parametric time-modulating function, in: 12th International Conference on Structural Safety and Reliability; Vienna, pp. 1–10.
  * Brune (1970) Brune, J.N., 1970.  Tectonic stress and the spectra of seismic shear waves from earthquakes.  Journal of geophysical research 75, 4997–5009.
  * Bruns (2004) Bruns, A., 2004.  Fourier-, hilbert-and wavelet-based signal analysis: are they really different approaches?  Journal of neuroscience methods 137, 321–332.
  * Cartwright and Longuet-Higgins (1956) Cartwright, D.E., Longuet-Higgins, M.S., 1956\.  The statistical distribution of the maxima of a random function.  Proceedings of the royal society of london. series a. mathematical and physical sciences 237, 212–232.
  * Castro and Munguia (1993) Castro, R., Munguia, L., 1993\.  Attenuation of p and s waves in the oaxaca, mexico, subduction zone.  Physics of the earth and planetary interiors 76, 179–187.
  * Chesneau (2023) Chesneau, C., 2023.  Extensions of two bivariate strict archimedean copulas.  Computational Journal of Mathematical and Statistical Sciences 2, 159–180.
  * Christen and Fox (2010) Christen, J.A., Fox, C., 2010\.  A general purpose sampling algorithm for continuous distributions (the t-walk).  Bayesian Analysis 5, 263–281.
  * Davy and Doucet (2003) Davy, M., Doucet, A., 2003.  Copulas: a new insight into positive time-frequency distributions.  IEEE signal processing letters 10, 215–218.
  * Ekström et al. (2012) Ekström, G., Nettles, M., Dziewoński, A., 2012.  The global cmt project 2004–2010: Centroid-moment tensors for 13,017 earthquakes.  Physics of the Earth and Planetary Interiors 200, 1–9.
  * Fraser and Boashash (1994) Fraser, G., Boashash, B., 1994\.  Multiple window spectrogram and time-frequency distributions, in: Proceedings of ICASSP’94. IEEE International Conference on Acoustics, Speech and Signal Processing, IEEE. pp. IV–293.
  * García et al. (2009) García, D., Singh, S.K., Herraiz, M., Ordaz, M., Pacheco, J.F., Cruz-Jiménez, H., 2009\.  Influence of subduction zone structure on coastal and inland attenuation in mexico.  Geophysical Journal International 179, 215–230.
  * Huang et al. (2022) Huang, Z., Xu, Y.L., Tao, T., 2022.  Multi-taper s-transform method for evolutionary spectrum estimation.  Mechanical Systems and Signal Processing 168, 108667.
  * Kramer (1996) Kramer, S.L., 1996.  Geotechnical earthquake engineering.  Pearson Education India.
  * Liang et al. (2007) Liang, J., Chaudhuri, S.R., Shinozuka, M., 2007.  Simulation of nonstationary stochastic processes by spectral representation.  Journal of Engineering Mechanics 133, 616–627.
  * Nelsen (2007) Nelsen, R.B., 2007.  An introduction to copulas.  Springer Science & Business Media.
  * Ordaz et al. (1994) Ordaz, M., Singh, S., Arciniega, A., 1994.  Bayesian attenuation regressions: an application to mexico city.  Geophysical Journal International 117, 335–344.
  * Parajuli (2017) Parajuli, R.R., 2017.  Earthquake magnitude, ground motion and structural safety (en inglés).  A Journal on Rural Infrastructure Development 8, 87–90.
  * Priestley (1965) Priestley, M.B., 1965.  Evolutionary spectra and non-stationary processes.  Journal of the Royal Statistical Society: Series B (Methodological) 27, 204–229.
  * Rosen et al. (2009) Rosen, O., Stoffer, D.S., Wood, S., 2009.  Local spectral analysis via a bayesian mixture of smoothing splines.  Journal of the American Statistical Association 104, 249–262.
  * Sánchez-Sesma et al. (2011) Sánchez-Sesma, F.J., Weaver, R.L., Kawase, H., Matsushima, S., Luzón, F., Campillo, M., 2011\.  Energy partitions among elastic waves for dynamic surface loads in a semi-infinite solid.  Bulletin of the Seismological Society of America 101, 1704–1709.
  * Sato (2013) Sato, T., 2013.  Fractal characteristics of phase spectrum of earthquake motion.  Journal of Earthquake and Tsunami 7, 1350010.
  * Shearer (2009) Shearer, P.M., 2009.  Introduction to seismology.  Cambridge university press.
  * Singh et al. (1989) Singh, S., Ordaz, M., Anderson, J., Rodriguez, M., Quaas, R., Mena, E., Ottaviani, M., Almora, D., 1989\.  Analysis of near-source strong-motion recordings along the mexican subduction zone.  Bulletin of the Seismological Society of America 79, 1697–1717.
  * Stafford et al. (2009) Stafford, P., Sgobba, S., Marano, G., 2009.  An energy-based envelope function for the stochastic simulation of earthquake accelerograms.  Soil Dynamics and Earthquake Engineering 29, 1123–1133.
  * Tang et al. (2016) Tang, Y., Lam, N., Lumantarna, E., Tsang, H.H., 2016\.  Generation of synthetic earthquake accelerograms based on up-to-date seismological ground motion models, in: Australian Earthquake Engineering Society 2016 Conference, Melbourne, pp. 1–14.
  * Vlachos et al. (2016) Vlachos, C., Papakonstantinou, K.G., Deodatis, G., 2016.  A multi-modal analytical non-stationary spectral model for characterization and stochastic simulation of earthquake ground motions.  Soil Dynamics and Earthquake Engineering 80, 177–191.
  * Vlachos et al. (2018) Vlachos, C., Papakonstantinou, K.G., Deodatis, G., 2018.  Predictive model for site specific simulation of ground motions based on earthquake scenarios.  Earthquake Engineering & Structural Dynamics 47, 195–218.
  * Wei (2006) Wei, W.W., 2006.  Time series analysis.