\documentclass[Journal,letterpaper, SectionNumbers]{ascelike-new}
%% Please choose the appropriate document class option:
% "Journal" produces double-spaced manuscripts for ASCE journals.
% "NewProceedings" produces single-spaced manuscripts for ASCE conference proceedings.
% "Proceedings" produces older-style single-spaced manuscripts for ASCE conference proceedings. 
%
%% For more details and options, please see the notes in the ascelike-new.cls file.

% Some useful packages...
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{lmodern}
\usepackage{graphicx}
\usepackage[figurename=Fig.,labelfont=bf,labelsep=period]{caption}
% \usepackage{subcaption}
\usepackage{amsmath}
%\usepackage{amsfonts}
%\usepackage{amsbsy}
\usepackage{mathtools}
\usepackage{amssymb}
\usepackage{mathrsfs}
\usepackage{newtxtext,newtxmath}
\usepackage[colorlinks=true,citecolor=red,linkcolor=black]{hyperref}
\usepackage{subfig}
\usepackage{multirow}
\usepackage{braket}
\mathchardef\mhyphen="2D
%

% Please add the first author's last name here for the footer:
\NameTag{Hexiang Wang, \today}
% Note that this is not displayed if the NoPageNumbers option is used
% in the documentclass declaration.
%
\begin{document}

% You will need to make the title all-caps
\title{Time  Domain  Intrusive  Probabilistic Seismic Risk Analysis of Nonlinear Shear Frame Structure}

\author[1]{Hexiang Wang}
\author[2]{Fangbo Wang}
\author[1]{Han Yang}
\author[1]{Yuan Feng}
\author[1]{Jeff Bayless}
\author[1]{Norman A. Abrahamson}
\author[1,3]{Boris Jeremi{\'c}*}


\affil[1]{Department of Civil and Environmental Engineering, University of California, Davis, CA, USA}
\affil[2]{School of Civil Engineering, Tianjin University, Tianjin, China}
\affil[3]{Earth and Environmental  Sciences Area,  Lawrence  Berkeley  National Laboratory, Berkeley, CA, USA, Email:jeremic@ucdavis.edu}

\maketitle

% Please include an abstract:
\begin{abstract}
Presented  is  a  time domain intrusive framework for probabilistic seismic risk
analysis.
% 
Seismic source characterization is mathematically formulated.
%
Methodology for simulating non-stationary seismic motions for given source, path
and site is proposed.
%
Both  uncertain motions and uncertain structural parameters are characterized as
random process/field and represented with Hermite polynomial chaos.
% 
Intrusive  modeling  of  Armstrong-Fredrick kinematic hardening based on Hermite
polynomial  chaos  is  formulated  and  incorporated  into  Galerkin  stochastic
elastic-plastic FEM.
%  
Time-evolving  probabilistic  structural  response  is  solved through developed
stochastic elastic-plastic FEM.
%
Following that, formulation for seismic risk analysis is derived.

The  framework  is  illustrated by seismic risk analysis of an eight-story shear
frame structure.
%
Uncertainties  are  propagated  from earthquake source into uncertain structural
system.
%
Difficulties  of  choosing  intensity  measure in the conventional framework are
avoided  since  all  the  uncertainties  and  important  characteristics  (e.g.,
spectrum  acceleration  $Sa$  and  peak  ground  acceleration  $PGA$) of seismic
motions are directly carried by the random process excitations in time domain.
%
Stochastic  dynamic  equations  are  solved  in  an intrusive way, circumventing
non-intrusive Monte Carlo simulations.
%
\end{abstract}

\begin{KeyWords}
Probabilistic seismic risk analysis \enspace Stochastic ground motion \enspace Stochastic elastoplastic FEM \enspace Time domain \enspace Fourier amplitude spectra
\end{KeyWords}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
Performance-based  Earthquake  Engineering  (PBEE)  \cite{Allin2000} has been a
successful  framework that allows for objective and quantitative decision-making
through seismic risk analyses.
%
% These decisions may be related to choosing between alternative designs, retrofitting options, or rehabilitation strategies, and the decision may be arrived at following economic, social or environmental considerations, or some combination of these. 
%
% Frameworks for PBEE, such as that of the Pacific Earthquake Engineering Research (PEER) center \cite{Allin2000}, allow this decision-making process to be followed in an objective and quantitative manner. 
%
Equation   \ref{eq:seismic_risk_analysis}   demonstrates   state   of   the  art
methodology of seismic risk analysis:

\begin{equation}
\lambda(EDP>z) = \int \underbrace{|\frac{d\lambda(IM>x)}{dx}|}_\text{PSHA} \underbrace{G(EDP>z|IM=x)}_{\text{fragility}} dx 
\label{eq:seismic_risk_analysis}
\end{equation}

\noindent 
where  $\lambda(EDP>z)$ is the annual rate of engineering demand parameter (EDP,
i.e., performance target) exceeding specific level $z$.
%
EDP  hazard is computed as convolution of probabilistic seismic hazard  analysis
(PSHA) results and structural fragility with respect to intensity
measure (IM) of ground shaking.
% 
PSHA,  usually  done  by  engineering seismologist, estimates exceedance rate of
intensity  measure $\lambda(IM>x)$ considering all possible faults and scenarios
near the engineering site.
%
Structural  fragility  $G(EDP>z|IM=x)$  defines the exceeding probability of EDP
given ground motion with particular IM level $x$.
%
With  properly defined damage measure (DM) as a function of EDP(s), seismic risk
of damage state can be calculated.

%
The  choice  of IM is crucial in seismic risk analysis, as it serves as proxy of
damaging  ground  motions and all the uncertainties in ground motion are assumed
could be represented by the variability of IM.
%
Spectral   acceleration  $Sa(T_0)$  is  commonly  adopted  as  IM  for  building
structures.
%
Many  ground  motion predictions equations (GMPEs) are developed to quantify the
median and aleatory variability of $Sa(T_0)$ \cite{Gregor2014}.
%
However,  the  problem  is  that  the  scalar spectral acceleration cannot fully
describe the influence of ground-motion variability upon engineering objects.
% 
\citeN{Stafford2010}  investigated  different  intensity measures and found that
they  are  generally  not strongly correlated, which indicates that knowledge of
just  one  IM  distribution  is  not  sufficient  to  describe  any of the other
ground-motion characteristics.
%


In addition, $Sa(T_0)$ as IM for surface building structures, is based on frequency domain, linear dynamic analysis of single degree of freedom system.
%
% Spectral acceleration, put forward based on dynamic response of single degree of freedom (SDOF) oscillator in frequency domain, is a very useful damaging proxy for linear elastic response of surface buildings, while it does not correlate well with time domain nonlinear response and are not suitable at all for any other types of engineering objects.    
% 
% It might not be appropriate for widespread nonlinear problems. 
%
When  nonlinear  inelastic and/or higher mode  response  is  expected,  use  of  $Sa(T_0)$  is  not
appropriate.
%
Nonlinear response history analysis (RHA) with spectrum-matched ground motion is
found   to   give   un-conservatively   biased  estimates  \cite{Iervolino2010,Huang2009}.
%
\citeN{Grigoriu2016}  showed  that  generally $Sa(T_0)$ is weakly dependent with
engineering  demand  parameters for realistic structures and fragilities defined
as functions of $Sa(T_0)$ have large uncertainties and of limited practical use.
%
% The bias phenomenon is attributed to asymmetric relationship between conditional spectral acceleration ordinates at periods affecting inelastic response \cite{Seifried2016}. 
%
Furthermore,  for  many  other  engineering objects (e.g., dams, deeply embedded
structures,  etc.),  it  is  very  difficult  to find a proper IM in engineering
practices.
%
For  example,  choice  of  IM  among peak ground acceleration (PGA), peak ground
velocity  (PGV), Arias intensity (AI) and cumulative absolute velocity (CAV) has
been   contentiously   argued   for   deformation  analysis  of  dam  embankment
\cite{Davoodi2013}.
%
Though  Vector-valued  PSHA  \cite{Baker2007}  was put forward to mitigate this
issue, it is rarely performed in practices.
%
The difficulty lies in fragility computation.
%
The  fragility  becomes  a function of vector IMs (e.g., a fragility surface for
two IMs), which requires a large number of structural analyses to be quantified.
%
Properly choosing multiple IMs is also a problem.
%
Many  times,  even if proper IMs, such as AI and CAV, are identified, additional
efforts are still needed to develop GMPE for these IMs and their correlation.
%

An  effective  solution  to  the  aforementioned  problems  would  be  to remove
intensity measure (IM) as an intermediate proxy from risk calculation.
%
% To this end, we established a time domain stochastic framework for seismic risk analysis.  
%
With  this  in mind, a time domain intrusive framework for probabilistic seismic
risk analysis is developed and described here.
%
The  framework  is  based  on  the  progress of Fourier amplitude spectrum (FAS)
modeling   of  seismic  motions  over  last  several  decades  \cite{Brune1970,Boore1983,Boore2003,Boore2015}.
%
Recent  advances  in  inter-frequency  correlation  of FAS \cite{Stafford2017b,Bayless2017}  and  Fourier phase derivative modeling \cite{Baglio2017} are also
taken into account.
%
Uncertain  motions  are simulated from stochastic FAS and Fourier phase spectrum
(FPS), and are modeled as non-stationary random process in time domain.
%
With   the   proposed  framework,  engineering  seismologists  do  not  need  to
interpret/simplify ground motion into IM(s).
%
Correspondingly,  structural  engineers  do  not need to compute fragility curve
based on IM.
%
Instead,  all the important characteristics and uncertainties in seismic motions
are   captured   through  the  random  process  and  propagated  into  uncertain
engineering   system   with   direct   ``communication''   between   engineering
seismologists and structural engineers.
%

Another  feature  of  the proposed framework is the circumvention of Monte Carlo
(MC) simulation.
%
MC  approach  is  non-intrusive  in  the  sense  that  no  modifications  to the
underlying deterministic solver are required.
%
The  state  of  probabilistic  space  is  characterized  by large, statistically
significant number of deterministic samplings of system random parameters.
%
In  conventional  seismic risk analysis, structural fragility curve is developed
by incremental dynamic analysis (IDA) \cite{Vamvatsikos2002}.
%
IDA,  though  theoretically straightforward, is numerically demanding because of
the slow convergence rate that is inherent in MC approach.
%
Hundreds  of  structural  analysis  need  to  be  performed  with deterministic
sampling  of  uncertain  material properties and uncertain ground excitations at
different IM levels.
%
The  same  issue  of  MC  approach  also limits the application of physics-based
seismic   waveform   modeling  techniques  \cite{Graves2010,Maechling2007}  into
hazard/risk analysis.
%
Millions  of  MC earthquake scenarios over regional geology have to be simulated
using    deterministic   wave   propagation   programs,   such   as   CyberShake
\cite{Graves2011}  considering uncertain kinematic sources, crustal geology and
site conditions.
%
\cite{Maechling2007}  estimated  that  ``it would require 300 million CPU-hours
and  well  over  100 years to complete all the simulations needed to calculate a
PSHA hazard curve''.
%
% Even    with    reciprocity    technique    developed    by    \citeN{Zhao2006},
% \citeN{Maechling2007}  states  that  ``it is still prohibitively computationally
% expensive to produce a waveform-based PSHA hazard map.''
%



To avoid non-intrusive MC simulation, Galerkin stochastic elastic-plastic finite element method
(SEPFEM) has been developed within the authors' research group over the years 
\cite{Jeremic2005a,Sett2005a,Sett2010a,Karapiperis2016,Wang2016,Wang2019}.     
%
Galerkin  SEPFEM  is  an  intrusive  approach,  requiring new developments based on
variational  formulation  of  the  underlying  stochastic  partial  differential
equations (SPDE).
%
Using  appropriate  choice  of  orthogonal  polynomial  chaos  basis,  intrusive
Galerkin  SEPFEM  guarantees optimal convergence rates, and is more efficient than
non-intrusive MC approach \cite{Xiu2010,Elman2011}.
%
Both  random field structural parameters and random  process seismic motions are
represented   by   Hermite   polynomial  chaos  (PC)  \cite{Sakamoto2002}  with
correlation structure characterized by Karhunen-Lo{\`e}ve (KL) expansion
\cite{Zheng2017}.
%
Using  Galerkin  SEPFEM, probabilistic  dynamic response of uncertain structural
system  driven  by  uncertain  seismic  motions  is  represented  by  unknown PC
coefficients.
%
Deterministic    linear  system  equations of these  unknown temporal-spatial PC
coefficients,  equivalent  to  the  original  stochastic  PDE,  are derived from
Galerkin projection technique in weak sense.
% 
Seismic risk is then computed from probabilistic dynamic structural response.
%

The organization of this paper is as follows: 
%
The  proposed  time  domain  intrusive  framework for probabilistic seismic risk
analysis is formulated in section \ref{sec:general_description}.
%
Next,  the proposed methodology is illustrated by a numerical example in section
\ref{sec:illustrative_example}     with    conclusions    drawn    in    section
\ref{sec:conclusions}.
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Time Domain Intrusive Framework for Seismic Risk Analysis}
\label{sec:general_description}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The  proposed  framework  consists  of  four  components,  as  shown  in  Figure
\ref{figure_risk_framework}:  seismic  source characterization (SSC), stochastic
ground  motion  modeling,  stochastic  finite  element analysis and seismic risk
computation.
%

\begin{figure}[!htbp]
\centering
\includegraphics[width=0.75\textwidth]{Figures/risk_framework_illustration.pdf}
\caption{\label{figure_risk_framework}   Time  domain  intrusive  framework  for
seismic risk analysis.}
\end{figure}

In  the  first  step,  SSC quantifies the uncertainty in earthquake scenarios so
that  the  probabilistic  scenario  space  $\lambda(M,  R,  \Theta)$ for a given
engineering  site  can  be  discretized  into  $N$  mutually exclusive events as
follows:
%

\begin{equation}
\lambda(M, R, \Theta) = \bigcup_{i=1}^{N} \lambda_i(M_i, R_i, \Theta_i)
\label{eq_discretize_scenario}
\end{equation}

\noindent 
where  $\lambda(\cdot)$  is the annual occurrence rate, $M$ is the magnitude and
$R$  is  the  distance metric, which could be either rupture distance $R_{rup}$,
hypocenter distance $R_{hyp}$, or Joyner-Boore distance $R_{jb}$.
%
$\Theta$  denotes  any  other  scenario metrics that are required for stochastic
ground  motion  modeling,  for example, style of fault, hanging wall identifier,
etc.
%
$N$  is  the total number of seismic scenarios considering all the active faults
in the region.
%
Basic  relations  for  seismic source characterization are formulated in Section
\ref{sec:SSC}.
%

For     each    scenario    event    $S_i(M_i,    R_i,    \Theta_i)$,    section
\ref{sec:stochastic_ground_motion_modeling}  presents  the procedure to simulate
time  domain uncertain motions from stochastic FAS and FPS using inverse Fourier
transform.
%
The   simulated   ground  motion  population  for  event  $S_i$  is  denoted  as
$\{\boldsymbol \Gamma_i\}$.
%

At  the  third  step, both uncertain motions and uncertain structural parameters
are   represented   by   Hermite   PC-KL  expansion  as  formulated  in  section
\ref{sec:pc_kl_expansion}.
%
Two  choices  are  provided  here: (1)  Random process characterization
(i.e.,  PC-KL  expansion)  is performed for  each individual motion population $\{\boldsymbol
\Gamma_i\}$  and  conduct  further  Galerkin  stochastic  FEM  analysis for each
scenario $S_i$.
%
(2) Seismic motion  population from  different  scenarios is first combined
as an ensemble population $\{\boldsymbol \Gamma\}$ following Equation \ref{eq_combine_population}:

\begin{equation}
\{ \boldsymbol \Gamma \} = \bigcup_{i=1}^{N} \{ w_i \otimes \boldsymbol \Gamma_i \}
\label{eq_combine_population}
\end{equation}

\noindent 
with 

\begin{equation}
w_i = \frac{\lambda_i}{\sum_{i=1}^{N} \lambda_i}
\label{eq_weight_population}
\end{equation}

\noindent 
where $\bigcup_{i=1}^{N} \{w_i \otimes \boldsymbol \Gamma_i \}$ denotes the weighted combination of population
$\{\boldsymbol   \Gamma_i   \}$   with   weight   $w_i$   defined   as  Equation
\ref{eq_weight_population}.
%
The  annual  occurrence rate of the ensemble population $\{\boldsymbol \Gamma\}$
is $\overline \lambda = \sum_{i=1}^{N} \lambda_i$.
%
The weighted combination can be performed by aggregating individual population $\{\boldsymbol{\Gamma}_i\}$ of different size $n_i$, $i=1, 2, ..., N$ such that $n_i$ is proportional to $w_i$, i.e., $w_i=n_i/\sum_{i=1}^{N}n_i$.
%
Clearly, size $n_i$ for all $i=1, 2, ..., N$ should be large enough to represent the random process motions from individual seismic scenario.
%
As a result, weighted ensemble population $\{\boldsymbol{\Gamma}\}$ with occurrence rate $\overline{\lambda}$ is statistically equivalent to the aggregation of motion population $\{\boldsymbol{\Gamma}_i\}$ from individual scenario with rate $\lambda_i$.
%
Then the ensemble population $\{\boldsymbol  \Gamma\}$ can be characterized as a
single   random   process   and  single stochastic  FEM  analysis  is  performed  with
PC-represented random process motions.
%
Compared with PC-KL representation for each individual population $\{\boldsymbol{\Gamma}_i\}$, the consequence of PC-KL expansion for ensemble population $\{\boldsymbol{\Gamma}\}$ is that larger dimension of PC is required since underlying random process of population $\{\boldsymbol{\Gamma}\}$ is more uncertain and less correlated among different times. 
%
If both individual population $\{\boldsymbol{\Gamma}_i\}$ and ensemble population $\{\boldsymbol{\Gamma}\}$ are accurately characterized by PC-KL expansion and propagated into uncertain structure through SFEM, EDP hazard can be calculated by either Equation 5 or Equation 6:
%
% Time-evolving  probabilistic structural response is obtained from stochastic FEM analysis.
%
% for  chosen  EDP(s),  failure  probability  $P_i(EDP>z|  \boldsymbol  \Gamma_i)$
% conditioned  on  individual  scenario  $S_i(M_i,  R_i,  \Theta_i)$  or $P(EDP>z|
% \boldsymbol \Gamma)$ conditioned on combined scenarios can be determined.
% %
% Finally, EDP hazard can be calculated     by     Equation
% \ref{eq:time_domain_risk_individual}     or    \ref{eq:time_domain_risk_combine}
% without using any IM(s) and performing Monte Carlo fragility simulations.

\begin{equation}
\lambda(EDP>z) = \sum_{i=1}^{N}\lambda_i(M_i, R_i, \Theta_i)P(EDP>z|\boldsymbol \Gamma_i)
\label{eq:time_domain_risk_individual}
\end{equation}
%
\begin{equation}
\lambda(EDP>z) = \overline \lambda P(EDP>z|\boldsymbol \Gamma)
\label{eq:time_domain_risk_combine}
\end{equation}
%
where $P_i(EDP>z|  \boldsymbol  \Gamma_i)$ is the failure probability conditioned on individual population $\{\boldsymbol{\Gamma}_i\}$ and $P(EDP>z| \boldsymbol \Gamma)$ is the failure probability conditioned on ensemble population $\{\boldsymbol{\Gamma}\}$.  
%
Both Equation 5 and Equation 6 give consistent result for EDP hazard. 
%
The difference is that by using Equation 5, many more less expensive SFEM analyses are performed while using Equation 6 requires a single, yet more expensive SFEM analysis. 
%
When the number of scenarios $N$ is small, it is practical to perform stochastic
FEM   analysis   for   each  scenario  and  compute  EDP hazard by  Equation
\ref{eq:time_domain_risk_individual}.
%
The  advantage  is  that  controlling  scenario  can  be identified through EDP hazard 
de-aggregation.
%
However, when there are many seismic scenarios, quantifying ensemble population as a single random process through PC-KL expansion and performing single stochastic FEM analysis can be computationally more efficient.
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
\subsection{Seismic Source Characterization}
\label{sec:SSC}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Seismic  source characterization (SSC) and earthquake rupture forecast (ERF) are
complex scientific issues.
%
Earthquake  occurrence  rate  tends  to be comprehensively evaluated by multiple
approaches,  for  example,  using  historical seismicity, geological information
(e.g.,  long term slip rates and paleoseismic recurrence intervals) and geodetic
information \cite{Field2017}.
%
Assuming  Poisson  process  of  earthquake  occurrence,  annual  occurrence rate
$\lambda^{f}$ of earthquakes on a fault can be estimated based on seismic moment
balance \cite{McGuire2004}:
%

\begin{equation}
\small
\lambda^{f} = \frac{\mu A \mathcal{S}}{\int_{0}^{M_{max}}E(M)f(M) \, dM}
\label{eq:activity_rate_moment_balancing1}
\end{equation}
%
\noindent 
where $\mathcal{S}$ is annual slip rate, $\mu$ is shear modulus of crust and $A$
is fault area, $f(M)$  is  the  probabilistic  model  of magnitude distribution, which could be
truncated  exponential  model,  Young's  and  Coppersmith  characteristic  model
\cite{Youngs1985}, truncated Gaussian model, etc.
%
The seismic moment of earthquake, $E(M)$ with magnitude $M$ is given as:
%

\begin{equation}
\small
E(M) = 10^{1.5M+16.05}
\label{eq:activity_rate_moment_balancing2}
\end{equation}

In  engineering  practices,  only  earthquakes  greater  than  certain magnitude
$M_{min}$ are considered, whose annual occurrence rate $\overline{\lambda}^{f} $
is:

%
\begin{equation}
\overline \lambda^{f} = \lambda^{f} \int_{M_{min}}^{M_{max}}f(M) dM
\end{equation}

%
Using  probabilistic  models  of rupture area conditioned on magnitude $f(A|M)$,
rupture  width conditioned on rupture area $f(W|A)$ \cite{Leonard2010}, rupture
location  along strike (AS) $f(Y)$ and down-dip (DD) $f(Z)$, distance metric $R$
and  other  scenario  metrics $\Theta$, for example, depth to the top of rupture
plane  $Z_{tor}$,  can be geometrically characterized as $g(R, \Theta| M)$ for a
given engineering site \cite{Hale2018}.
%
The  discretized mutually exclusive scenarios $\lambda_i(M_i, R_i, \Theta_i)$ in
Equation \ref{eq_discretize_scenario} is then quantified as:
%

\begin{equation}
\lambda_i(M_i, R_i, \Theta_i) = \sum_{j=1}^{m} \overline{\lambda_j}^{f} \int_{\boldsymbol{\Lambda_i}} f_j(M) g_j(R, \Theta|M) \,dM\,dR\,d\Theta
\end{equation}

\noindent 
where  $m$  is  the  total  number  of  active faults, subscript $j$ denotes the
probabilistic   models   and   quantities   specific   to  the  $j^{th}$  fault,
$\boldsymbol{\Lambda_i}$ is the integral domain for the $i^{th}$ discretized
scenario  with  magnitude  step $\Delta M$, distance step $\Delta R$ and $\Delta
\Theta$ for any other scenario metrics $\Theta$ if required:
%
\vspace{-0.4cm}

\begin{equation}
\boldsymbol \Lambda_i = [M_i-\frac{\Delta M}{2}, M_i + \frac{\Delta M}{2}] \times [R_i-\frac{\Delta R}{2}, R_i + \frac{\Delta R}{2}] \times [\Theta_i - \frac{\Delta \Theta}{2}, \Theta_i + \frac{\Delta \Theta}{2}]  
\label{eq:intergal_domain}
\end{equation}

%
Many PSHA programs could perform SSC, e.g., HAZ45 \cite{Hale2018}.
%
It is noted that presented above are fundamental relations for seismic characterization of fault sources.
%
% The formulations for area sources are similar and omitted.
%
For regions with unknown fault locations or having background seismicity, areal source should also be considered and characterized. 
%
See references \cite{Coppersmith2012,Moschetti2015} for more details on seismic source characterization of areal source.   
%
Epistemic  uncertainties  in  slip rate, magnitude distribution models and other
parameters,   which   are   typically   considered   with  logic  tree  approach
\cite{Musson2012}, are not considered here for simplicity.
%
In  addition, for some sites, authoritative estimates of magnitude, location and
rate  of  earthquake  ruptures  could  be  determined  from established regional
earthquake  rupture forecast (ERF) models, for example, UCERF3 \cite{Field2017}
for California region.
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
\subsection{Time Domain Stochastic Ground Motion Modeling}
\label{sec:stochastic_ground_motion_modeling}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stochastic  modeling of FAS and FPS are two indispensable components to simulate
% time domain uncertain motions.
%
Time domain uncertain motions can be simulated from stochastic FAS and Fourier phase derivative \cite{Boore2003b,Boore2003}. 
%
Specifically,  uncertain  FAS  of  seismic  motions  is  modeled  as  Log-normal
distributed  random  field  \cite{Bora2015,Stafford2017b}  in  frequency space,
whose  marginal  median  behavior  is  simulated  by  the  stochastic  method of
\citeN{Boore2003}.
%
It is referred to as Boore03 approach hereafter.
%
Boore03   approach   simulates   FAS   using   $w^2$  radiated  source  spectrum
\cite{Brune1970}  with  modification  for  path  and  site effects, as shown in
Equation \ref{eq:FAS_SIMSM}:
%

\begin{equation}
FAS(f) = A_0(M_0, f)Z(R)exp(-\pi f R/Q\beta)S(f)exp(-\pi \kappa_0 f)
\label{eq:FAS_SIMSM}
\end{equation}

%
\noindent 
where $M_0$ is the seismic moment; $\beta$ is the source shear wave velocity;
%
$Z(R)$  and $exp(-\pi f R/Q\beta)$ represent the contribution from path effects:
$Z(R)$ is the geometrical spreading term as a function of distance $R$.
%
Term $exp(-\pi  f  R/Q\beta)$  quantifies the anelastic attenuation as the inverse of
the regional quality factor, $Q$.
%
The  site effects including site amplification through crustal velocity gradient
and  near  surface  attenuation are demonstrated by $S(f)$ and
$\kappa_0$ filter $exp(-\pi \kappa_0 f)$, respectively.
%
Term  $A_0$    represents   the  radiated  acceleration  source  spectrum, which
could be characterized by single-corner-frequency model:
% ) or double-corner-frequency model \cite{Boore2014b}.  
%

\begin{equation}
A_0(M_0, f) = CM_0\bigg[\frac{(2\pi f)^2}{1+(f/f_0)^2}\bigg]
\label{eq:single_corner_frequency_model}
\end{equation}

\noindent 
where $f_0$ is the corner frequency, which in Brune's model \cite{Brune1970} is
related to source stress drop $\Delta \sigma$ as follows:
%
% $C$ is a scaling constant, typically taken as $R_{\Theta \phi}VF/(4\pi\rho\beta^3)$, in which $R_{\Theta \phi}$ is the average radiation coefficient 0.55 for shear waves \cite{Boore1984}, $V=1/\sqrt{2}$ is the partitioning of seismic energy into horizontal components.  
%
% $F=2$ is factor for free surface reflection.  

\begin{equation}
f_0 = 4.9\times10^6\beta (\Delta \sigma/M_0)^{1/3}
\end{equation}

Boore03  approach  is  well-recognized  for  its simplicity and effectiveness to
capture the marginal mean behavior of stochastic FAS.
%
% It   is   noted   that   several   GMPEs  of  FAS  have  been  derived  recently
% \cite{Bora2015,Bora2018,Bayless2018b}  and  are  also  suitable  to  model  the
% marginal behavior of FAS within the proposed framework.
%
\citeN{Bayless2018c}  pointed out that the inter-frequency correlation structure
of FAS random field is also important for seismic risk analysis.
%
Misrepresenting-representing   the   correlation   structure,   e.g.,   assuming
inter-frequency independence, would lead to underestimation of seismic risk.
%
% Recently, inter-frequency correlation models for $\epsilon$, which is defined as residual of uncertain $FAS(f)$ with respect to its median estimate normalized by standard derivation, have been developed \cite{Stafford2017b,Bayless2017}. 
%
Therefore,  inter-frequency  correlation  model  for  stochastic  FAS  developed
recently \cite{Stafford2017b,Bayless2017} is adopted here.
%
% Adopted here is the FAS correlation structure identified by \cite{Bayless2017}.
%

Though  the behavior  of  FAS was well studied, modeling Fourier phase
angles is still challenging.
%
% In  the  original  Boore03  approach  \cite{Boore2003},  random  phase  info is
% contained  in  the  stationary  Gaussian  white  noise  modulated by an envelope
% function.
%
Conventionally random  phase  info is simulated using stationary  Gaussian  white  noise  modulated by an envelope function.
%
However,  \citeN{Montaldo2003}  stated  that  conventional  Gaussian white noise approach could not reliably reproduce the non-stationarity of ground motions.
%
For  this  reason,  the  use of phase difference $\Delta \Phi$  was suggested by \citeN{Ohsaki1979}.
%
Using California strong ground motion data, \citeN{Thrainsson2002} modeled phase
differences as Beta distribution.
%
However,  the  established  phase  difference  models are affected by the signal
length of each record.
%
It  is  more stable to normalize phase difference by signal length and study the
probabilistic    model    of    phase   derivative   $\dot{\Phi}$   defined   as
\cite{Boore2003b} :
%

\begin{equation}
\dot{\Phi} = \frac{\Delta \Phi}{\Delta f}
\label{eq:phase_derivative_model}
\end{equation}

%
Based   on   3551   ground   motion  records  from  PEER  NGA-West  1  database,
\citeN{Baglio2017}   found   that   the  distribution  of  phase  derivative  is
leptokurtic and fits well to Logistic model:
%

%%%% new equation 
\begin{equation}
f(\dot{\Phi}; \mu, \sigma) = \frac{1}{4\sigma} sech^{2}(\frac{\dot{\Phi}-u}{2\sigma})
\label{eq:phase_derivative_logistic_model}
\end{equation}

\noindent 
where  $\mu$  and  $\sigma$  are  the  mean  and scale parameter of the Logistic
distribution  $f(\dot{\Phi};  \mu,  \sigma)$,  $sech(\cdot)$  is  the hyperbolic
secant function.
%
Following  \citeN{Baglio2017},  the  mean  value  $\mu$  is a fixed parameter to
position the distribution along the signal length.
%
For example, setting mean parameter $\mu$ equal to $\pi/df$ would align the peak of uncertain seismic motions to the center of simulated signal length. 
%
The  prediction equation of scale parameter $\sigma$ is correlated to earthquake
magnitude  $M$,  rupture  distance  $R_{rup}$,  $V_{s30}$  and directivity index
$D_{Dir}=R_{hyp}  - R_{rup}$ with coefficients $\alpha_1$, $\alpha_2$, $\beta_1$
$\sim$  $\beta_4$,  $\gamma_1$ and $\gamma_2$ determined from maximum likelihood
estimation:

\begin{equation}
log(\sigma/\pi) = \alpha_1 + \alpha_2 log[\beta_1 + 10^{\beta_2M} + \beta_3 R_{rup} + \beta_4 log(V_{s30})+ \gamma_1 + \gamma_2 D_{Dir}]
\label{eq:GMPE_phase_derivative}
\end{equation}

Phase  derivatives  $\dot{\Phi}(f)$  among  frequency  coordinates is modeled as
Logistic  distributed  random  field following exponential correlation 
with correlation length $l_f=0.05Hz$:

%
\begin{equation}
\displaystyle
{Cov(\dot{\Phi}(f_1), \dot{\Phi}(f_2)) = e^{-\displaystyle{\frac{|f_1-f_2|}{l_f}}}}
\label{eq:correlation_structure_phase_derivative}
\end{equation}

% \begin{equation}
% \displaystyle
% {Cov(\dot{\Phi}(f_1), \dot{\Phi}(f_2)) = e^{-\frac{|f_1-f_2|}{l_f}}}
% \label{eq:correlation_structure_phase_derivative}
% \end{equation}
%
The  methodology  of time domain stochastic ground motion modeling is summarized
below:
%
\begin{enumerate}

\item  Compute  marginal  median of Log-normal distributed random field $FAS(f)$
following Boore03 approach.
%


\item  Generate  realizations  of  Log-normal  distributed random field $FAS(f)$
according  to  the marginal estimation in step 1 and inter-frequency correlation
model by \citeN{Bayless2017}.
%


\item  Determine  the  scale  parameter  $\sigma$ of marginal Logistic model for
phase derivative random field with Equation \ref{eq:GMPE_phase_derivative}.
%
Set mean value $\mu$ to $\pi/df$ for central peak \cite{Baglio2017}.
%


\item Generate realizations of Logistic distributed random field $\dot{\Phi}(f)$
with marginal distribution from step 3 and exponential correlation structure.
%


\item  Multiply  realization of phase derivative in step 4 by frequency interval
$df$ to get realizations of phase difference $\Delta \Phi(f)$.
%
Compute  realizations  of  phase angles from $\Delta \Phi(f)$. First phase angle
can be randomly set and would not affect the synthesis.
%


\item Inverse Fourier transformation to time domain with realizations of FAS and
FPS.


\end{enumerate}

Time domain stochastic ground motion for a single scenario with magnitude $M=7$,
distance $R_{rup}=15km$ has been simulated following the above methodology.
%
The  detailed  modeling  parameters  for  source,  path and site effects of this
scenario               are              given              in              Table
\ref{table:stochastic_ground_motion_modeling_parameter}.
%
The   marginal   median  FAS  is  computed  using  program  SMSIM  developed  by
\citeN{Boore2005}.
%
% Finite faults effects are considered through equivalent point source distance $R_{PS}$ following \cite{Boore2015}. 
%
With       reference      to      recent      GMPE      studies      of      FAS
\cite{Bora2015,Bora2018}, marginal lognormal standard deviation of
FAS  has been adopted as total $\sigma=0.8 \,ln$ units.
%
The maximum modeling frequency $f_{max}$ is 20Hz.
%
It is noted that ergodic assumption was used in developing these  GMPEs of FAS.
%
A smaller value of marginal standard deviation can be used for non-ergodic probabilistic seismic risk analysis if additional source, path or site specific information is available. 
%

\begin{table}[!htbp]
\footnotesize
\centering
% \caption{\label{table:stochastic_ground_motion_modeling_parameter} Input parameters for median FAS modeling with SMSIM \cite{Boore2005}}
\caption{\label{table:stochastic_ground_motion_modeling_parameter}  Source, path
and  site  parameters  for stochastic ground motion modeling of seismic scenario
$M=7$, $R_{rup}=15km$.}
\begin{tabular}{ccc}
\hline
\multicolumn{1}{l}{Parameter type} & Name                    & Value                                                                                                                                               \\ \hline
\multirow{6}{*}{Source}            & Magnitude               & M=7                                                                                                                                                 \\
                                   & Source density          & $\rho_s=2.8g/cm^3$                                                                                                                                  \\
                                   & Source velocity         & $\beta =3.6km/s$                                                                                                                                    \\
                                   & $w^2$ source spectrum   & Single corner frequency with $\Delta \sigma=8.0 MPa$                                                                                              \\
                                   & Fault type              & Reverse fault $F_{rv}=1$                                                                                                                            \\
                                   & Dip angle               & $45^{\circ}$                                                                                                                                        \\ \hline
\multirow{4}{*}{Path}              & Distance metrics        & \begin{tabular}[c]{@{}c@{}}$R_{rup}=15km$, $R_{hyp}=18km$\\ $R_{jb}=12km$, $R_x=-12km$\end{tabular}                                                 \\
                                   & Finite faults effects   & \begin{tabular}[c]{@{}c@{}}Equivalent point source model \cite{Boore2015} \\ with $R_{PS}=22.18km$\end{tabular}                                                                           \\
                                   & Geometrical spreading   & \begin{tabular}[c]{@{}c@{}}Hinged line segments model \cite{Atkinson1995} \end{tabular} \\
                                   & Anelastic attenuation Q & Three line segments model by  \cite{Boore2003}                                                                                                                      \\ \hline
\multirow{2}{*}{Site}              & Site amplification      & \begin{tabular}[c]{@{}c@{}}$V_{s30}=620m/s$\\  Table 4 of \cite{Boore2015} \end{tabular}                                                        \\
                                   & $\kappa_0$ attenuation    & $\kappa_0=0.03s$                                                                                                                                    \\ \hline
\end{tabular}
\end{table}

Combining  stochastic  FAS  with  uncertain  Fourier  phase  info,  500 realizations of time  domain stochastic motions are synthesized.
%
Figure \ref{figure_acc_realizations} shows three different synthesized   accelerations. 
%
Large  variability is observed, for example, peak ground acceleration could vary
from 1.8$m/s^2$ to 5.5$m/s^2$.

\begin{figure}[!htbp]
\centering
\includegraphics[width=0.31\textwidth]{Figures/Acc_time_series_realiztion70.pdf}
\includegraphics[width=0.29\textwidth]{Figures/Acc_time_series_realiztion100.pdf}
% \includegraphics[width=0.31\textwidth]{Figures/Acc_time_series_realiztion350.pdf}
\includegraphics[width=0.325\textwidth]{Figures/Acc_time_series_realiztion367.pdf}
\caption{\label{figure_acc_realizations}  Realizations of uncertain acceleration
time series population.}
\end{figure}

Spectral  acceleration ($Sa$) of 500 synthesized realizations are calculated and
compared   with   weighted   average   prediction   of  five  NGA  West-2  GMPEs
\cite{Gregor2014}  with weights 0.22 for ASK14, 0.22 for BSSA14, 0.22 for CB14,
0.22 for CY14 and 0.12 for I14.
%
From  figure  \ref{figure_SA_verification_GMPE}(a), median spectral acceleration
$Sa$ from simulated stochastic ground motion is in very good agreement with GMPE
predictions for all period ordinates.
%
No systematic bias is observed.
%
Figure  \ref{figure_SA_verification_GMPE}(b) shows that the  standard  deviation  of simulated response spectra is around 0.65 $ln$ units, which is
consistent with aleatory variability of $Sa$ given by GMPEs.
%
In   other   words,   time  domain  stochastic  ground  motions  simulated  with
aforementioned  methodology  could  not only characterize the median behavior of
$Sa$  very well, but also carry desired amount of uncertainty that is consistent
with empirical GMPEs.

\begin{figure}[!htbp]
\centering
\subfloat[GMPE verification of $Sa$]{
\includegraphics[width=0.51\textwidth]{Figures/SA_GMPE_verification_std_08_no_smooth.pdf}}
\subfloat[Bias check]{
\hspace{-0.3cm}
\includegraphics[width=0.5\textwidth]{Figures/Standard_deviation_std_08_no_smooth.pdf}}
\caption{\label{figure_SA_verification_GMPE}     Verification    of    simulated
stochastic  motions: (a) Check median spectral acceleration $Sa$ with NGA West-2 GMPE
(b) Check aleatory variability of simulation spectral acceleration $Sa$.}
\end{figure}
%
The marginal distribution of simulated accelerations at all the time instances is observed to be Gaussian.
%
Similar observation is also made by \citeN{Wang2016} from statistical analysis of seismic records.
%
Therefore, time domain stochastic ground motions is modeled as a Gaussian distributed non-stationary random process.
%
The  random  process  would  be  represented  with  Hermite  polynomial chaos as
formulated in the next section.
%

It  is  worthwhile  to  mention  that  the random process incorporates much more
information about uncertain ground motions than GMPE used in conventional PBEE.
%
GMPE  only  quantifies  the  variability of selected IM, such as $Sa$, while the
random  process  carries  not  only  consistent variability of $Sa$ but also any
other important characteristics, e.g., PGA, CAV, etc.
%
Realistic inter-frequency correlation of FAS \cite{Bayless2017} is captured.
%
Non-stationarity of seismic motions is quantified through phase derivative modeling without using any modulation function.
%
Compared  with  existing  ground  motion modeling techniques commonly adopted by
reliability  community, e.g., evolutionary power spectrum and white noise random
phase    spectrum,   presented   methodology   is   directly   compatible   with
state-of-the-art seismic source characterization.
%
It  could  explicitly  account for specific source, path and site condition in
both stochastic modeling of FAS and FPS.
%
Many  reliability  analysis  methods,  such  as  probabilistic density evolution
method \cite{Xu2019} can be readily combined with the presented methodology and
incorporated into the proposed risk analysis framework for PBEE.
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Hermite   Polynomial   Chaos Karhunen-Lo{\`e}ve expansion}
\label{sec:pc_kl_expansion}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
This section formulates Hermite polynomial chaos Karhunen-Lo{\`e}ve (PC-KL)
expansion for general heterogeneous random field $D(\boldsymbol{x}, \theta)$ of
arbitrary marginal distribution.
%
Both uncertain motions and uncertain structure parameters can be represented
with PC-KL expansion.
%
To achieve this, we first represent heterogeneous random field
$D(\boldsymbol{x}, \theta)$ of arbitrary marginal distribution through Hermite
polynomial chaos of underlying Gaussian heterogeneous random field
$\gamma(\boldsymbol{x}, \theta)$ up to order $P$ \cite{Sakamoto2002,Wang2016}:
%
% The  populations  of  simulated  time  series  seismic  motions  are regarded as
% realizations of underlying non-stationary random process.
%
% This random process can be quantified using Polynomial Chaos Karhunen-Lo{\`e}ve (PC-KL) expansion.
%
% Any heterogeneous random field  of arbitrary marginal
% distribution (Gaussian or non-Gaussian) can  be  represented  through  Hermite polynomial chaos (PC) up to
% order $P$ as \cite{Sakamoto2002,Wang2016}:
%

\begin{equation}
D(\boldsymbol{x}, \theta) 
= 
\sum_{i=0}^{P} D_i(\boldsymbol{x})\Omega_i(\gamma(\boldsymbol{x}, \theta)) 
\label{eq:PC_expansion}
\end{equation} 
% 

\noindent 
% where  $\boldsymbol{x}$  is  the  general  coordinate in
% space $V$ that can span space and/or time.
%
where $\theta$ denotes the  uncertainties.
%
Functions  $\{\Omega_i\}=  \{1, \gamma, \gamma^2-1, \gamma^3-3\gamma, ...\}$ are
orthogonal,  zero  mean ($i\geq 1$) Hermite PC bases constructed from zero  mean,  unit variance kernel Gaussian
random field $\gamma(\boldsymbol{x}, \theta)$.
%
Then at the second step, Gaussian random field $\gamma(\boldsymbol{x}, \theta)$ can be further decomposed by Karhunen-Lo{\`e}ve
(KL) theorem \cite{Zheng2017}. 

%
The  deterministic  PC coefficient field $D_i(\boldsymbol{x})$ can be calculated
through   marginal  distribution  of  $D(\boldsymbol{x},  \theta)$,   as shown
in Equation
\ref{eq:PC_coefficient}, where $\braket{\cdot}$ is the expectation operator.
%
\begin{equation}
D_i = \frac{\braket{D\Omega_i}}{\braket{\Omega_i^2}}
\label{eq:PC_coefficient}
\end{equation}
%
% The main idea behind Equation \ref{eq:PC_expansion} is to represent arbitrary heterogeneous random field of any type of marginal distribution with a kernel Gaussian random field.
% %
% In this way, deterministic PC coefficients field can be determined and 
%
The  covariance  structure  of  the  original  random field $Cov_D(x_1, x_2)$ is
mapped to the Gaussian covariance kernel $Cov_{\gamma}(x_1, x_2)$ as:
%
\begin{equation}
Cov_D(x_1, x_2) = \sum_{i=1}^{P} D_i(x_1) \; D_i(x_2)\;  i\,! \; Cov_{\gamma}(x_1, x_2) 
\label{eq:correlation_transformation}                        \; 
\end{equation}
%
The  Gaussian covariance kernel $Cov_{\gamma}(x_1, x_2)$ can be eigen-decomposed
into  probabilistic spaces up to dimension $M$, according to Karhunen-Lo{\`e}ve
(KL) theorem \cite{Zheng2017}:

%
\begin{equation}
\gamma(\boldsymbol{x}, \theta) 
= 
\sum_{i=1}^{M} \sqrt{\lambda_i} f_i(\boldsymbol{x}) \xi_i(\theta)
\label{eq:KL_expansion}
\end{equation}
%
%
where  $\{\xi_i(\theta)\}$  are  the multidimensional, orthogonal, zero mean and
unit variance Gaussian random variables, and
$\lambda_i$  and $f_i(\boldsymbol{x})$ are the eigen-values and eigen-vectors of
the  covariance  kernel  $Cov_{\gamma}(x_1,  x_2)$ that satisfy Fredholm's integral
equation of the second kind \cite{Sakamoto2002}.
%

% \begin{equation}
% \int_V Cov_{\gamma}(x_1, x_2) f_i(x_1) \, dx_1 = \lambda_{i} f_i(x_2)
% \end{equation}
%
% Due  to  the  unit variance of Gaussian random field $\gamma(\boldsymbol{x},
% \theta)$,  eigen-values $\lambda_i$.  and eigen-vectors $f_i(\boldsymbol{x})$,
% of
% the  covariance  kernel  $Cov_{\gamma}(x_1,  x_2)$, 
% are constrained as follows:
% %
% \begin{equation}
% \sum_{i=1} (\sqrt{\lambda_i} f_i(\boldsymbol{x}))^2 = 1
% \label{eq:constrained_KL_condition}
% \end{equation}
% %



Combining  Equations \ref{eq:PC_expansion}  and  \ref{eq:KL_expansion}, the
resultant  PC-KL  representation  of  general  random  field  $D(\boldsymbol{x},
\theta)$  is obtained as,
%
%
\begin{equation}
D(\boldsymbol{x}, \theta) = \sum_{i=0}^{K} d_i(\boldsymbol{x}) \Psi_i(\{\xi_j(\theta)\})
\label{eq:PC-KL expansion}
\end{equation}
%
%
%  \noindent 
%  {\huge \bf NOTE: 
%  where  $\{\Psi_i\}$   are multi-dimensional orthogonal Hermite PC bases of order
%  $P$, dimension $M$, probabilistic space (i.e., \{$\xi_j(\theta)$\}).}
%  %
%  %HEXIANG upper sentence needs to be rewritten, last part is not clear
%
%===== Rephrased as =====
where  $\{\Psi_i\}$   are multi-dimensional orthogonal Hermite PC bases of order
$P$ constructed from $M$ dimensional probabilistic space (i.e., \{$\xi_j(\theta)$\}, $j=1, 2, ..., M$).
%
%
The total number of multidimensional Hermite PC bases $K$ is related to order $P$
and dimension $M$ as $K = 1+\sum_{s=1}^{P}\frac{1}{s!}\prod_{j=0}^{s-1}(M+j)$. 

%
% \begin{equation}
% \displaystyle K = 1+\sum_{s=1}^{P}\frac{1}{s!}\prod_{j=0}^{s-1}(M+j)
% \end{equation}
%

By  equating  two  representations  of  $D(\boldsymbol{x}, \theta)$ in Equations
\ref{eq:PC_expansion}   and   \ref{eq:PC-KL   expansion},  the  coefficients  of
multi-dimensional Hermite PC are derived as:
%

% \begin{equation}   
% d_i(\boldsymbol{x}) = \frac{P!}{\braket{\Psi_i^2}} D_P(\boldsymbol{x}) \prod_{j=1}^{P} \frac{\sqrt{\lambda_{k(j)}}f_{k(j)}(\boldsymbol{x})}{\sqrt{\sum_{m=1}^{M}(\sqrt{\lambda_m}f_m(\boldsymbol{x}))^2}}
% \end{equation}
%
\begin{equation}
d_i(\boldsymbol{x}) = \frac{p!}{\braket{\Psi_i^2}} D_p(\boldsymbol{x}) \prod_{j=1}^{p} \frac{\sqrt{\lambda_{k(j)}}f_{k(j)}(\boldsymbol{x})}{\sqrt{\sum_{m=1}^{M}(\sqrt{\lambda_m}f_m(\boldsymbol{x}))^2}}
\end{equation}
%
where $p$ is the order of the polynomial $\Psi_i$. 
%
From  Equation  \ref{eq:PC-KL  expansion},  PC  synthesized  marginal  mean  and variance of the original heterogeneous random field can be calculated as:
%
\begin{equation}
\braket{D(\boldsymbol{x}, \theta)} = d_0(\boldsymbol{x})
\label{eq:PC_marginal_mean}
\end{equation}
%
\begin{equation}
Var(D(\boldsymbol{x}, \theta)) = \sum_{i=1}^{K} d_i^2(\boldsymbol{x}) \braket{\Psi_i^2}
\label{eq:PC_marginal_variance}
\end{equation}
%
PC-synthesized correlation structure can also be computed as:

\begin{equation}
\displaystyle
Cov_D(x_1, x_2) = \frac{\sum_{i=1}^{K} d_i(x_1)d_i(x_2)\braket{\Psi_i^2}}{\sqrt{Var(D(x_1))Var(D(x_2))}}
\label{eq:PC_correlation_structure}
\end{equation}
%

Equations    \ref{eq:PC_marginal_mean},      \ref{eq:PC_marginal_variance}   and
\ref{eq:PC_correlation_structure}  can be used to compare the PC-synthesized
statistics  with statistics of original random field $D(\boldsymbol{x}, \theta)$
and check the goodness of PC-KL expansion.
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Galerkin Stochastic Finite Element Method}\label{sec:SFEM}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Stochastic Galerkin approach intrusively solves the stochastic PDE with optimal convergence \cite{Ghanem1996,Sett2010a,Wang2016}.
%
% This is different from the non-intrusive methods, such as Monte Carlo and stochastic collocation \cite{Babuska2007b}. 
%
Stochastic  Galerkin approach intrusively solves the stochastic partial
differential equations (PDE) with optimal convergence
\cite{Sett2010a,Wang2016}.
%
Compared  to  deterministic  finite element method (FEM),  Galerkin  stochastic 
FEM introduces spectral discretization of probabilistic domain in addition to
the spatial and temporal discretization.
%
Using  standard  spatial  FEM  discretization, unknown displacement random field
$u(\boldsymbol{x},   t,   \theta)$   can   be   expressed  with  shape  function
$N_i(\boldsymbol{x})$ and uncertain displacement $u_i(t, \theta)$ at nodes:

%
\begin{equation}
u(\boldsymbol{x}, t, \theta) = \sum_{i=1}^{N} N_i(\boldsymbol{x}) u_i(t, \theta)
\label{eq:spatial_discretization}
\end{equation}
%
Uncertain  displacement at nodes $u_i(t, \theta)$, can be further represented with
aforementioned  multidimensional Hermite PC basis $\phi_j(\{\xi_r(\theta)\})$ of
dimension $M^{u}$, order $P^{u}$:
%

\begin{equation}
u_i(t, \theta) = \sum_{j=0}^{K^{u}} u_{ij}(t) \phi_j(\{\xi_r(\theta)\})
\label{eq:probabilistic_discretization}
\end{equation}
%
Combining          Equations         \ref{eq:spatial_discretization}         and
\ref{eq:probabilistic_discretization},     spatial-probabilistic     discretized
expression of $u(\boldsymbol{x}, t, \theta)$ is given:
%
\begin{equation}
u(\boldsymbol{x}, t, \theta) = \sum_{i=1}^{N}\sum_{j=0}^{K^{u}} N_i(\boldsymbol{x})  u_{ij}(t) \phi_j(\{\xi_r(\theta)\})
\end{equation}
%
Galerkin weak formulation of stochastic partial differential equations of motion
can then be written in the following form:
%
\begin{center}
\begin{eqnarray}
\sum_e \bigg[ 
& & 
\int_{D_e} N_m(\boldsymbol x) \rho(\boldsymbol x) N_n(\boldsymbol x) d\,V \, \ddot{u}_n(t, \theta) 
+ 
\nonumber
\\
&+&  
\int_{D_e} B_m(\boldsymbol{x}) E(\boldsymbol{x, \theta}) B_n(\boldsymbol x) d\,V\,u_n(t, \theta) - f_m(t, \theta) \bigg] 
=
0
% \end{split}
\label{eq:weak_form_FEM}
\end{eqnarray}
\end{center}
%
where  $\displaystyle  \sum_e$  denotes  the  assembly procedure over all finite
elements, while $\rho(\boldsymbol{x})$ is deterministic material density field.
%
The shape function gradient function $B_n(\boldsymbol x)$ is given as:

%
\begin{equation}
B_n(\boldsymbol x) = \nabla N_n(\boldsymbol x)
\end{equation}
%
In Equation~\ref{eq:weak_form_FEM}, 
$E(\boldsymbol{x}, \theta)$ is uncertain tangential stiffness matrix,  while
$f_m(t, \theta)$ is uncertain global force vector that incorporates various
elemental contributions.
%
% This uncertain force vector can be used to input uncertain seismic excitations,
% using Domain Reduction Method \cite{Bielak2001,Yoshimura2001}.
% %

Expansion of  uncertain stiffness matrix $E(\boldsymbol{x},  \theta)$,  and
uncertain force vector $f_m(t, \theta)$ into Hermite PC bases
$\Psi_k(\{\xi_r(\theta)\})$   and   $\psi_l(\{\xi_r(\theta)\})$   of   dimension
$M^{E}$, order $P^{E}$ and dimension $M^{f}$, order $P^{f}$, respectively, yields:
%
%
\begin{equation}
E(\boldsymbol{x}, \theta) 
= 
\sum_{k=0}^{K^{E}} E_k(\boldsymbol{x}) \Psi_k(\{\xi_r(\theta)\})
\label{eq:stiffness_discretization}
\end{equation}
%
%
\begin{equation}
f_m(t, \theta)= \sum_{l=0}^{K^{f}} f_{ml}(t) \psi_l(\{\xi_r(\theta)\})
\label{eq:force_discretization}
\end{equation}
%
By combining equations            \ref{eq:probabilistic_discretization},
\ref{eq:stiffness_discretization}    and    \ref{eq:force_discretization} and 
equation \ref{eq:weak_form_FEM}, one obtains:
%

%
\begin{equation}
\begin{split}
& \sum_e \bigg[ \int_{D_e} N_m(\boldsymbol x) \rho(\boldsymbol x) N_n(\boldsymbol x) d\,V \, \sum_{j=0}^{K^u} \ddot{u}_{nj} \phi_j(\{\xi_r(\theta)\} - \sum_{l=0}^{K^{f}} f_{ml} \psi_l(\{\xi_r(\theta)\}) \\
& + \int_{D_e} B_m(\boldsymbol{x}) \sum_{k=0}^{K^{E}} E_k(\boldsymbol{x}) \Psi_k(\{\xi_r(\theta)\}) B_n(\boldsymbol x) d\,V\, \sum_{j=0}^{K^u} u_{nj} \phi_j(\{\xi_r(\theta)\}\bigg] = 0
\end{split}
\label{eq:weak_formulatio_expand}
\end{equation}
%

By performing  Galerkin  projection of Equation \ref{eq:weak_formulatio_expand}
onto  PC  bases  $\phi_i(\{\xi_r(\theta)\})$, to minimize the residual, system of
deterministic   ordinary   differential   equations   (ODE)  involving  temporal
derivative  of  unknown  PC  coefficients  $u_{nj}$,  is  developed:
% with Einstein notation:
%

\begin{equation}
M_{minj} \ddot{u}_{nj} + K_{minj} u_{nj} = F_{mi}
\label{eq:deterministic_ODE}
\end{equation}
%
where mass tensor/matrix $M_{minj}$ is given by equation~\ref{eq:Mass_stochastic_matrix}:
%

%  {\huge \bf NOTE: 
%  Hexiang, since material density is determinstic, this mass matrix is then
%  deterministic! Is this right?}
%  %

%%%%%% === Hexiang Reply: Yes, because of the determinisitc material density, the expanded mass matrix is only determinisitc FEM mass matrix multiplied by the double product of PC bases. And because of the orthogonality of PC bases, it ends up with diagonally duplicate of deterministic mass matrix. 

\begin{equation}
M_{minj} = \sum_{e} \int_{D_e} N_m(\boldsymbol x) \rho(\boldsymbol x) N_n(\boldsymbol x) d\,V \, \braket{\phi_i \phi_j} 
\label{eq:Mass_stochastic_matrix}
\end{equation}
%
stochastic stiffness tensor/matrix $K_{minj}$  is given by equation~\ref{eq:stiffness_stochastic_matrix}:
%
\begin{equation}
K_{minj} = \sum_{k=0}^{K^{E}} \sum_{e} \int_{D_e} B_m(\boldsymbol{x})  E_k(\boldsymbol{x})  B_n(\boldsymbol x) d\,V\, \braket{\Psi_k \phi_i \phi_j}
\label{eq:stiffness_stochastic_matrix}
\end{equation}
%
and stochastic force tensor/vector  $F_{mi}$ is given by
equation~\ref{eq:force_stochastic_matrix}

\begin{equation}
F_{mi} = \sum_{l=0}^{K^{f}} f_{ml} \braket{\psi_l \phi_i} 
\label{eq:force_stochastic_matrix}
\end{equation}


In                   Equations                  \ref{eq:Mass_stochastic_matrix},
\ref{eq:stiffness_stochastic_matrix}    and    \ref{eq:force_stochastic_matrix},
terms $\braket{\phi_i        \phi_j}$,        $\braket{\psi_l       \phi_i}$       and
$\braket{\Psi_k\phi_i\phi_j}$ are the ensemble average tensors of double-product
and tri-product of different PC bases. 
%
These  ensemble  average tensors could be pre-computed and used to construct the
stochastic mass matrix $M_{minj}$ and stochastic stiffness matrix $K_{minj}$.
%
It is noted that Einstein's notation for
tensor indices summation is assumed throughout \cite{local-63}.
%

The deterministic system of ordinary differential equations (ODE) from Equation
\ref{eq:deterministic_ODE}, can be integrated in time using 
dynamic integrator algorithms, for example Newmark method \cite{Newmark1959},
or Hilber-Hughes-Taylor $\alpha$-method \cite{Hilber1977}.
%
Result of such time marching solution will be time histories of displacement PC
coefficients $u_{nj}$.
%
Those time evolving  displacement  PC coefficients $u_{nj}$ can then be used 
to develop complete probabilistic dynamic finite element response.
%
With resulting complete probabilistic dynamic finite element response, 
any damage  measure, in fact all damage measures related to  EDP(s) can be
applied  to  trace the failure probability $P_i(EDP>z|\boldsymbol{\Gamma}_i)$ or
$P(EDP>z|\boldsymbol{\Gamma})$.
%
EDP hazard   can    then    be    computed    according    to    Equations
\ref{eq:time_domain_risk_individual} and \ref{eq:time_domain_risk_combine}.



The  above formulation of Galerkin stochastic FEM is complete for linear elastic
problem  with  constant  uncertain  elastic  stiffness matrix $E(\boldsymbol{x},
\theta)$.
%
For  nonlinear,    inelastic    problems,  additional  formulation of stochastic
elastic-plastic  FEM  (SEPFEM) is required and relies on recent  developments
\cite{Jeremic2005a,Sett2005a,Sett2009b,Sett2010a,Arnst2012,Rosic2014,Karapiperis2016}.
% 
One   of  the  challenges  of  the    SEPFEM    lies    in    the development of
the  probabilistic  elastic-plastic stiffness at the constitutive level that is to be
used for finite element level computations.
%
Eulerian-Lagrangian form of the Fokker-Planck-Kolmogorov (FPK) equation has been
successfully    used    to    obtain       probabilistic      stress   solutions
\cite{Jeremic2005a,Sett2005a,Sett2010a}.
%  
It  is  noted  in order to produce uncertain stiffness, least       square      
optimization and linearizion techniques
\cite{Sett2010a,Karapiperis2016}   are used.
%


To             this      end,    in one dimension (1D), elastic plastic material
model    with    vanishing    elastic    region  is  used  in  conjunction  with
Armstrong-Fredrick nonlinear kinematic hardening \cite{Armstrong66,Dettmer2004}.
%
This  approach simplifies modeling, as elastic plastic response directly follows
Armstrong-Frederick nonlinear equation.
%
% For approach proposed here, direct modeling of nonlinear, probabilistic inelastic
% inter-story drift, the probabilistic nonlinear response between
% inter-story  restoring  force
% $F^{R}$  and  inter-story  drift  $\eta$  is    formulated through direct  PC-based 
% Galerkin intrusive 
% probabilistic     modeling    of    Armstrong-Fredrick    hysteretic    behavior.
%
For the approach proposed here, probabilistic nonlinear response between inter-story restoring force $F^{R}$ and inter-story drift $\eta$ is formulated through direct PC-based Galerkin intrusive probabilistic modeling of Armstrong-Fredrick hysteretic behavior.



In   incremental   form,   Armstrong-Fredrick   kinematic   hardening   relation
\cite{Armstrong66} between inter-story  restoring  force $F^{R}$ and 
inter-story  drift $\eta$ can
be written as:

\begin{equation}
dF^{R} = H_a\, d\eta - C_r F^{R} |d\eta|
\label{eq:AF_constitutive}
\end{equation}
%
where $H_a$ and $C_r$ are model parameters. 
%
By      setting    $dF^{R}    =0$,  the  ultimate  inter-story  restoring  force
becomes $F^{R}_{max} =H_a/C_r$.
%
The    tangential    stiffness    $E(F^{R})$   is a  function of restoring force
$F^{R}$:
%

\begin{equation}
E(F^{R}) = \frac{d F^{R}}{d \eta} = H_a - C_r F^{R}\,sgn(d\eta) 
\label{eq:tangential_stiffness}
\end{equation}
%
where $sgn(\cdot)$ is the sign function. Equation~\ref{eq:tangential_stiffness} 
can be written as:

\begin{equation}
E(F^{R}) =  H_a \pm C_r F^{R} 
\label{eq_simplified_tangential_stiffness}
\end{equation}
%
where   $+$ sign is taken for negative inter-story  drift $d\eta$ and $-$ sign is
taken for positive inter-story drift $d\eta$.
%
In  the  general  probabilistic  setting, model parameters $H_a$  and $C_r$ can be
uncertain and
modeled  as  random  fields $H_a(\boldsymbol x, \theta)$ and $C_r(\boldsymbol x,
\theta)$.
%
By applying  PC  expansion with Hermite PC bases $\varphi_i(\{ \xi_r(\theta)\})$ to
those two model parameters, the following equations are obtained:
%
%
\begin{equation}
H_a(\boldsymbol{x}, \theta) = \sum_{i=0}^{P} H_{ai}(\boldsymbol x) \varphi_i(\{ \xi_r(\theta)\})
\label{eq:Ha_expanded}
\end{equation}
%
%
\begin{equation}
C_r(\boldsymbol{x}, \theta) = \sum_{i=0}^{P} C_{ri}(\boldsymbol x) \varphi_i(\{ \xi_r(\theta)\})
\end{equation}
%
%
The    inter-story  drift  increments  $d   \eta(\boldsymbol{x},  \theta)$, that
represent input to the to constitutive driver (Equation
\ref{eq:AF_constitutive}) are also  uncertain due to the  probabilistic
structural response $u(\boldsymbol x, t, \theta)$:
%
%
\begin{equation}
d \eta(\boldsymbol{x}, \theta) = \sum_{i=0}^{P} d \eta_i(\boldsymbol x) \varphi_i(\{ \xi_r(\theta)\})
\end{equation}
%
As  a  result,  probabilistic incremental restoring force $dF^{R}(\boldsymbol x,
\theta)$  and probabilistic tangential stiffness $E(\boldsymbol{x}, \theta)$ are
then:
%
%
\begin{equation}
dF^{R}(\boldsymbol x, \theta) = \sum_{i=0}^{P} dF^{R}_i(\boldsymbol x) \varphi_i(\{ \xi_r(\theta)\})
\end{equation}
%
%
\begin{equation}
E(\boldsymbol{x}, \theta) = \sum_{i=0}^{P} E_i(\boldsymbol x) \varphi_i(\{ \xi_r(\theta)\})
\label{eq:E_expanded}
\end{equation}
%
%
Substituting   Equations   \ref{eq:Ha_expanded} $\sim$ \ref{eq:E_expanded} 
into              Equations             \ref{eq:AF_constitutive}             and
\ref{eq_simplified_tangential_stiffness}  and  applying Galerkin projection on
PC basis $\varphi_i\{ \xi_r(\theta)\}$ yields:
%
%
\begin{equation}
\sum_{m=0}^{P} dF_m^{R} \braket{\varphi_m \varphi_i} = \sum_{j=0}^{P} \sum_{k=0}^{P} H_{aj}d\eta_{k} \braket{\varphi_j \varphi_k \varphi_i} \pm  \sum_{l=0}^{P} \sum_{n=0}^{P} \sum_{s=0}^{P} C_{rl} F^{R}_{n} d\eta_{s} \braket{\varphi_l \varphi_n \varphi_s \varphi_i}
\end{equation}
%
%
\begin{equation}
\sum_{i=0}^{P} E_m \braket{\varphi_m \varphi_i} = \sum_{j=0}^{P} H_{aj} \braket{\varphi_j \varphi_i} \pm  \sum_{l=0}^{P} \sum_{n=0}^{P} C_{rl} F^{R}_{n} \braket{\varphi_l \varphi_n \varphi_i}
\end{equation}
%
%
% By using  the  orthogonality  of  Hermite  PC  bases $\braket{\varphi_i \varphi_j}=
% \delta_{ij}$,   solutions   to   the  unknown  PC  coefficients  of  incremental
% inter-story force $dF^{R}(\boldsymbol x,
% \theta)$    and  inter-story  stiffness  $E(\boldsymbol{x},   \theta)$  can  be 
% written as:
%
By using  the  orthogonality  of  Hermite  PC  bases $\braket{\varphi_i \varphi_j}=
0$ for $i\neq j$, solutions   to   the  unknown  PC  coefficients  of  incremental
inter-story force $dF^{R}(\boldsymbol x,
\theta)$    and  inter-story  stiffness  $E(\boldsymbol{x},   \theta)$  can  be 
written as:
%
% \begin{equation}
% dF_i^{R} = \frac{1}{\braket{\varphi_i^2}} \bigg[ H_{aj}d\eta_{k} \braket{\varphi_j \varphi_k \varphi_i} \pm C_{rl} F^{R}_{n} d\eta_{s} \braket{\varphi_l \varphi_n \varphi_s \varphi_i} \bigg]
% \end{equation}
%
\begin{equation}
dF_i^{R} = \frac{1}{\rm Var[\varphi_i]} \bigg[ H_{aj}d\eta_{k} \braket{\varphi_j \varphi_k \varphi_i} \pm C_{rl} F^{R}_{n} d\eta_{s} \braket{\varphi_l \varphi_n \varphi_s \varphi_i} \bigg]
\end{equation}
%
% \begin{equation}
% E_i = H_{ai} \pm \frac{1}{\braket{\varphi_i^2}} C_{rl} F^{R}_{n} \braket{\varphi_l \varphi_n \varphi_i}
% \end{equation}
%
\begin{equation}
E_i = H_{ai} \pm \frac{1}{\rm Var[\varphi_i]} C_{rl} F^{R}_{n} \braket{\varphi_l \varphi_n \varphi_i}
\end{equation}

%
where $\braket{\cdot}$ is the expectation operator.
%
$\rm Var[\varphi_i]$ is the scalar variance of PC basis $\varphi_i\{
\xi_r(\theta)\}$, which equals to $\braket{\varphi_i^2}$. 
%
It is noted that in the above equations, Einstein's tensor summation notation is used with index $i$ as a free index. 
%
Terms  $\braket{\varphi_j  \varphi_k  \varphi_i}$,  $\braket{\varphi_l \varphi_n
\varphi_i}$  and  $\braket{\varphi_l  \varphi_n  \varphi_s  \varphi_i}$  are the
expectation   of   triple   and  quadruple  product  of  PC  basis  $\varphi_i\{
\xi_r(\theta)\}$.
%

The  above  1D  formulation  for  SEPFEM  is  implemented  in  the  context  of 
explicit, forward Euler algorithm, 
%
The  expanded stiffness  matrix  $K_{minj}$  is constructed  using stiffness  PC
coefficients      $^{(n)}E_i$     at     step     $n$     following     Equation
\ref{eq:stiffness_stochastic_matrix}.
%
Displacement   PC  coefficients  $^{(n+1)}u_{nj}$  of step $n+1$ are then solved
by  applying force vector $^{(n+1)}F_{mi}$ and using stiffness matrix $K_{minj}$
within Equation \ref{eq:deterministic_ODE}.
%
Following  that, incremental inter-story drift PC coefficients $^{(n+1)}d\eta_i$
are  calculated  from  displacement  response  $^{(n+1)}u_{nj}$  and incremental
uncertain restoring force $^{(n+1)}dF_i^{R}$ can be quantified as:
%
% %
% \begin{equation}
% {}^{(n+1)}dF_i^{R} 
% = 
% \frac{1}{\braket{\varphi_i^2}} \bigg[ H_{aj}\,^{(n+1)}d\eta_{k} 
% \braket{\varphi_j \varphi_k \varphi_i} \pm C_{rl} \, {}^{(n)}F^{R}_{n} 
% \, 
% {}^{(n+1)}d\eta_{s} \braket{\varphi_l \varphi_n \varphi_s \varphi_i} \bigg]
% \end{equation}
% %
%
\begin{equation}
{}^{(n+1)}dF_i^{R} 
= 
\frac{1}{\rm Var[\varphi_i]} \bigg[ H_{aj}\,^{(n+1)}d\eta_{k} 
\braket{\varphi_j \varphi_k \varphi_i} \pm C_{rl} \, {}^{(n)}F^{R}_{n} 
\, 
{}^{(n+1)}d\eta_{s} \braket{\varphi_l \varphi_n \varphi_s \varphi_i} \bigg]
\end{equation}
%
Updating the restoring force $^{(n+1)}F_i^{R}$ is then: 

\begin{equation}
{}^{(n+1)}F_i^{R} = {}^{(n)}F_i^{R} + {}^{(n+1)}dF_i^{R} 
\end{equation}
%
%
while new stiffness PC coefficients $^{(n+1)}E_i$ at step $n+1$ are then:
%
%
% \begin{equation}
% {}^{(n+1)}E_i 
% = 
% H_{ai} \pm \frac{1}{\braket{\varphi_i^2}} C_{rl} 
% \, 
% {}^{(n+1)}F^{R}_{n} \braket{\varphi_l \varphi_n \varphi_i}
% \end{equation}
%
\begin{equation}
{}^{(n+1)}E_i 
= 
H_{ai} \pm \frac{1}{\rm Var[\varphi_i]} C_{rl} 
\, 
{}^{(n+1)}F^{R}_{n} \braket{\varphi_l \varphi_n \varphi_i}
\end{equation}


%
%
%Other  iterative  algorithms,  such  as Backward Euler \cite{Jeremic96b}, could
%also be formulated similarly as in deterministic elastic-plastic FEM and are not
%repeated for conciseness.
%







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Illustrative Example}
\label{sec:illustrative_example}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


To  illustrate  the  proposed  framework,  seismic risk of a typical eight story
shear   frame   structure that has been studied by many researchers \cite{Li2006b,Mitseas2018,Papazafeiropoulos2017,Xu2019} , is developed. 
%
% Structure used in this study has   been  studied  earlier by \citeN{Xu2019}.
%
The frame structure is shown in Figure \ref{figure_eight_story_building}.
% 
%
\begin{figure}[!htbp]
\centering
\includegraphics[width=10.5cm]{Figures/shear_frame_illustration_update.pdf}
\caption{\label{figure_eight_story_building} Eight-story shear frame structure with uncertain floor stiffness under non-stationary seismic motions.}
\end{figure}

The hysteretic restoring force versus inter-story drift behavior is described by
Armstrong-Fredrick model presented in section \ref{sec:SFEM}.
%
Material parameter   $H_a$   of  Armstrong-Fredrick  model  is  assumed  to  be 
Gaussian distributed random field with $15\%$ coefficient of variation.
%
Means  of material parameter $H_a$ are given for different floors as: 
$H_{a1} \sim H_{a2}$ $1.59\times 10^7 N/m$,
$H_{a3}  \sim H_{a6}$ $1.66\times 10^7 N/m$ and 
$H_{a7} \sim H_{a8}$ $1.76\times 10^7 N/m$.
%
The      correlation  structure  of parameter $H_a$ is assumed to be exponential
between different
floors, with correlation length of $l_c=10$ floors.
%
Material parameter $C_r$ is assumed to be $C_r = 17.6$~1/m.
%
The resultant mean hysteretic behavior of first floor is also shown in Figure \ref{figure_eight_story_building}. 
%
Floor masses are assumed to be deterministic.
%
Rayleigh  damping  $\boldsymbol C=\alpha \boldsymbol M + \beta \boldsymbol K$ is
used  with  parameters  $\alpha$  and  $\beta$  chosen to be $\alpha=0.22$Hz and
$\beta=0.008s$.
%
Other structure     modeling     parameters     are     given     in     Table
\ref{tab_structure_parameter}.
%
Those parameters are determined from \citeN{Xu2019}. 
%
Parameters $H_a$ and $C_r$ are calibrated to match the hysteretic behavior shown in \citeN{Xu2019}.   

\begin{table}[!htbp]
\caption{\label{tab_structure_parameter}Parameters   of  the  eight-story  shear
frame structure.}
\centering
\small
\begin{tabular}{cccccc}
%\hline
$h_0[m]$ & $h[m]$ & $m_1$$\sim$$m_2[kg]$ & $m_3$$\sim$$m_4[kg]$ & $m_5$$\sim$$m_6[kg]$  & $m_7$$\sim$$m_8[kg]$ \\ 
\hline
$4$    & $3$         & 2$\times 10^{5}$      & 2.2$\times 10^{5}$  & 2.4 $\times 10^{5}$  & 2.5 $\times 10^{5}$ \\ 
%\hline
\end{tabular}
\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Seismic Source Characterization}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


The structure is
located  at  coordinate  ($10km$,  $40km$), $50km$ away from a strike slip fault
with    $90^{\circ}$    dip    angle,    as        shown       in      Figure   
\ref{figure_single_fault_problem_setup}.
%
%
\begin{figure}[!htbp]
\centering
\includegraphics[width=0.60\textwidth]{Figures/single_fault_problem_setup.pdf}
\caption{\label{figure_single_fault_problem_setup}  Seismic  risk analysis of an
eight-story  shear  frame  structure (red triangle) with uncertain stiffness $K$
subjected to earthquakes from a strike slip fault (black line).}
\end{figure}
%
The fault length is $250km$ with annual slip rate of $20mm/yr$.
%
Detailed  geometry  and  model  parameters  for SSC of the strike slip fault are
given in Table \ref{table:parameter_single_fault}.
%
%
\begin{table}[!htbp]
\caption{\label{table:parameter_single_fault}   Parameters  for  seismic  source
characterization (SSC) of the strike slip fault.}
\centering
\footnotesize
\begin{tabular}{cc}
%\hline
Parameter                             & Value                                                                      \\ 
\hline
Fault length                          & 250km                                                                      \\
Fault width                           & 15km                                                                       \\
Dip angle                             & $90^{\circ}$                                                                         \\
Slip rate $\mathcal{S}$               & 20mm/yr                                                                    \\
Style of faulting                     & Strike slip                                                                \\
% Coordinate                          & (60,-100), (60, 150)                                                      \\
$f(M)$                                & Truncated normal with $\sigma$=0.2 $n\sigma_{max}$=2 \cite{Hale2018}\\
$f(A|M)$                              & Delta function at $log(A) = M-4$                                             \\
$f(W|A)$                              & Delta function at $W=\sqrt{1.5A}$, limited to fault width  \\
$f(Y)$                                & Uniform distribution                                         \\
$f(Z)$                                & Uniform distribution                                                       \\ 
%\hline
\end{tabular}
\end{table}
%
Mean   characteristic  magnitude  of  the  fault $\overline{M}$    is   7.6, and
is related to fault area $A$ \cite{Leonard2010} as:
%
\begin{equation}
\small
\overline{M} = log(A) + 4
\label{eq:mean_characteristic}
\end{equation}
%
%
Only  earthquakes  with  magnitude  greater  than 5 (i.e. $M_{min}=5$) are
considered.
%
Following  the  procedure  of  SSC  in  section  \ref{sec:SSC},  annual  rate of
earthquakes occurring on the fault is $\overline{\lambda}=0.0067$/yr.
%
Probabilistic  scenario  space  $\lambda(M, R, \theta)$ is discretized into four
mutually  exclusive  scenario events $S_i(M_i, R_i, \Theta_i)$ as shown in Table
\ref{table:seismic_scenarios}.
%
The  computation is performed with probabilistic seismic hazard analysis program
HAZ45  \cite{Hale2018}  using  0.2  for  magnitude  step $\Delta M$ and 2km for
distance step $\Delta R$.
%
\begin{table}[!htbp]
\caption{\label{table:seismic_scenarios}  Seismic  scenarios for the strike slip
fault.}
\small
\centering
\begin{tabular}{cccc}
\hline
\ \ Scenario  ID & $M$   & $R_{rup}$ [km]     & Annual rate $\lambda(M, R_{rup})$ \\ \hline
\ \ \ 1&\ \ \ \ 7.3\ \ \ \ & 56  & 9.54$\times 10^{-4}$ \\
\ \ \ 2&\ \ \ \ 7.5\ \ \ \ & 56    & 2.40$\times 10^{-3}$ \\
\ \ \ 3&\ \ \ \ 7.7\ \ \ \ & 56    & 2.40$\times 10^{-3}$ \\
\ \ \ 4&\ \ \ \ 7.9\ \ \ \ & 56    & 9.54$\times 10^{-4}$ \\ \hline
\end{tabular}
\end{table}
%
%
% Activity rate $N(M>M_{min})$ is computed as by balancing moment rate on the single fault.  
%
% Combining the above activity rate $N(M>M_{min})$ with probabilistic models of rupture dimension and hypocenter location, full seismic source characterization can be performed with many available hazard analysis programs. 
%
% A list of seismic scenarios with different magnitude $M$, distance $R$ and annual occurrence rate $\lambda(M, R)$ can be obtained. 
%
% In this study, program HAZ45 \cite{Hale2018} was used for SSC of the strike slip fault. 
%
% Developed are 4 different scenarios as shown in Table \ref{table:seismic_scenarios}. 
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Time Domain Stochastic Ground Motion Modeling and Representation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
For  each  characterized  seismic  scenario  $S_i(M_i,  R_i,  \Theta_i)$,  $500$
realizations  $\{\boldsymbol  \Gamma_i\}$  of  time domain uncertain motions are
simulated    using    methodology   described   in      section
\ref{sec:stochastic_ground_motion_modeling}.
%
Figure \ref{figure_time_series_colletction} shows the first 200 realizations of simulated motions for earthquake scenario 1 with $M=7.3$, $R_{rup}=56km$.  

\begin{figure}[!htbp]
\centering
\includegraphics[width=0.6\textwidth]{Figures/Acc_realization_200.pdf}
\caption{\label{figure_time_series_colletction} Realizations of uncertain seismic motions for scenario $M=7.3$, $R_{rup}=56km$.}
\end{figure}

In  this  study,  ground  motion  populations  from four different scenarios are
combined   into   a   single   population   $\boldsymbol{\Gamma}$  using Equation
\ref{eq_combine_population} and modeled as a non-stationary random process.
%
The  random process is represented by multi-dimensional Hermite polynomial chaos
(PC) following the technique formulated in section \ref{sec:pc_kl_expansion}.
%
Since  marginal  distribution  of  the random process is observed to be Gaussian
(section    \ref{sec:stochastic_ground_motion_modeling}),    theoretically,    PC
representation with order 1 is sufficient.
%
The  dimension  of  PC  basis  needs  to  be carefully chosen to reconstruct the
correlation structure of the original random process.
%
To ensure the accuracy of PC-KL representation, following error measurements are
defined and evaluated:
%

\begin{itemize}

\item  
The absolute error on marginal mean of the random process:

\begin{equation}
% \varepsilon_m = \frac{1}{N_t} \sum_{k=1}^{N_t}|\frac{\mu(t_k)-\hat{\mu}(t_k)}{\sigma(t_k)}|
\varepsilon_m = \frac{1}{N_t} \sum_{k=1}^{N_t}|\mu(t_k)-\hat{\mu}(t_k)|
\end{equation}
%


\item  
The absolute error on marginal standard deviation of the random process:

\begin{equation}
\varepsilon_{std} = \frac{1}{N_t} \sum_{k=1}^{N_t}|\sigma(t_k)-\hat{\sigma}(t_k)|
\end{equation}
%

\item  
The absolute error on correlation of the random process: 

\begin{equation}
\varepsilon_{corr} = \frac{1}{N_t^2} \sum_{k=1}^{N_t} \sum_{l=1}^{N_t}|Cov(t_k, t_l)-\hat{Cov}(t_k, t_l)|
\end{equation}
%

\end{itemize}


\noindent
where  $\mu(t_k)$,  $\sigma(t_k)$  and  $Cov(t_k,  t_l)$  are the marginal mean,
marginal  standard  deviation  and  correlation field of simulated ground motion
population $\boldsymbol \Gamma$.
%
Terms $\hat{\mu}(t_k)$,  $\hat{\sigma}(t_k)$  and $\hat{Cov}(t_k, t_l)$ are statistics
calculated   from   PC   representation  of  the  random  process from Equations
\ref{eq:PC_marginal_mean},           \ref{eq:PC_marginal_variance}           and
\ref{eq:PC_correlation_structure}.
%
Term $t_k$  denotes  the $k^{th}$ time instance and $N_t$ is the total number of
time instances.
%




Hermite  PC  bases  of  order  1, dimension 20, 70, 150 and 300 are examined for
PC-KL expansion of random process motions.
%
The      errors    for      different      PC      bases   are  given  in  Table
\ref{tab_error_measurements}.
%
%
%
\begin{table}[!htbp]
\footnotesize
\caption{\label{table_motion_characterization_error}   Error   in  probabilistic
characterization  of non-stationary acceleration and displacement random process
using PC-KL expansion with different dimensions.}
\centering
\label{tab_error_measurements}
\begin{tabular}{ccccc}
%\hline
Dimension of PC                & Dim. 20 & Dim. 70 & Dim. 150 & Dim. 300 \\ \hline
Displacement mean error $\varepsilon_m$       & 8.63$\times 10^{-9}$ & 8.63$\times 10^{-9}$ & 8.63$\times 10^{-9}$  & 8.63$\times 10^{-9}$ \\
Displacement S.D. error $\varepsilon_{std}$       & 1.28$\times 10^{-7}$ & 1.28$\times 10^{-7}$ & 1.28$\times 10^{-7}$  & 1.28$\times 10^{-7}$  \\
Displacement correlation error $\varepsilon_{corr}$ & 0.059   & 2.26$\times 10^{-4}$ & 8.27$\times 10^{-6}$  & 3.06$\times 10^{-7}$  \\ \hline
Acceleration mean error $\varepsilon_m$       & 9.84$\times 10^{-9}$ & 9.84$\times 10^{-9}$  & 9.84$\times 10^{-9}$   & 9.84$\times 10^{-9}$   \\
Acceleration S.D. error $\varepsilon_{std}$       & 1.23$\times 10^{-7}$ & 1.23$\times 10^{-7}$ & 1.23$\times 10^{-7}$  & 1.23$\times 10^{-7}$  \\
Acceleration correlation error $\varepsilon_{corr}$ & 0.185   & 0.091   & 0.053    & 0.028    \\ 
%\hline
\end{tabular}
\end{table}
%
%
It  can  be  observed that  in  all  the four cases marginal behavior of the random
process  motions  is well captured with very small magnitudes of errors $\varepsilon_m$ and
$\varepsilon_{std}$.
%
As shown in Figure \ref{figure_marginal_verification}, synthesized marginal mean
and  marginal  standard  deviation  from  PC representation match very well with
statistics of simulated motions.
%
% Due to the transient nature of uncertain seismic motions, the marginal mean of the random process is close to 0.  
%
\begin{figure}[!htbp]
\centering
\subfloat[Displacement]{
\hspace{-1.4cm}
\includegraphics[width=0.6\textwidth]{Figures/KL_dis_verification.pdf}}
\subfloat[Acceleration]{
\hspace{-0.8cm}
\includegraphics[width=0.6\textwidth]{Figures/KL_acc_verification.pdf}}
\caption{\label{figure_marginal_verification}  Comparison between PC-synthesized
(black  line)  marginal mean and marginal standard deviation (SD) and statistics
of simulated ground motion realizations (red line).}
\end{figure}





As  the  dimension  of PC increases, the relative error of correlation structure
decreases while the computational cost in stochastic FEM increases.
%
It  is  noted that PC dimension 70 is already adequate to capture the relatively
smooth random displacement correlation field.
%
However,  acceleration  correlation  field  synthesized  from PC dimension 70 is
overestimated among many time steps.
%
PC  dimension  150  and  300 approximate acceleration correlation structure much
better.
%
Eventually,  considering  both  accuracy  and efficiency, Hermite PC of order 1,
dimension  150 is used to spectrally discretize the random process seismic
motions in stochastic FEM analysis.
%
The  comparison  between  the exact correlation structure and the PC synthesized
correlation structure is shown in Figure \ref{figure_motion_correlation}.
%
%
\begin{figure}[!htbp]
\centering
\vspace{-5mm}
\subfloat[Exact displacement correlation field]{
\includegraphics[width=0.50\textwidth]{Figures/KL_corr_exact_dis.pdf}}
%\subfloat[PC synthesized displacement correlation field]{
\subfloat[Synthesized displacement correlation field]{
\includegraphics[width=0.50\textwidth]{Figures/KL_corr_PC_dis.pdf}}
\vspace{-5mm}
\\
%\vspace{-5mm}
\subfloat[Exact acceleration correlation field]{
\includegraphics[width=0.50\textwidth]{Figures/KL_corr_exact_acc.pdf}}
%\subfloat[PC synthesized acceleration correlation field]{
\subfloat[Synthesized acceleration correlation field]{
\includegraphics[width=0.50\textwidth]{Figures/KL_corr_PC_acc.pdf}}
\vspace{-2mm}
\caption{\label{figure_motion_correlation}     Verification    of   PC synthesized
acceleration and displacement correlation field with PC dimension 150.}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Stochastic Galerkin FEM Analysis and Seismic Risk}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In order to perform stochastic Galerkin   FEM  analysis, it is also necessary to
characterize the randomness of stiffness of the structural system.
%
In  order  to  do  that, Hermite  PCs  of  dimension  2,  4 and 6 is used for
capturing  the  exponential  correlation  structure  of  random  field parameter
$H_a(\boldsymbol x, \theta)$.
%
%
%
%
It can be observed from  Figure \ref{figure_KL_stiffness} that PC dimension 4 can
reasonably well reconstruct the correlation of $H_a(\boldsymbol x, \theta)$.
%
\begin{figure}[!htbp]
\centering
\includegraphics[width=0.5\textwidth]{Figures/KL_verification_stiffness.pdf}
\caption{\label{figure_KL_stiffness} Characterization of exponential correlation
(correlation   length   $l_c=10$   floors)  of  uncertain  structural  parameter
$H_a(\boldsymbol x, \theta)$ using PCs of different dimensions.}
\end{figure}
%
% As mentioned before, parameter $C_r$ is assumed to be deterministic. 
%
% With  PC characterized $H_a$, the probabilistic hysteretic behavior of restoring
% force  versus  inter-story  drift  can  be  intrusively  modeled  following  the
% stochastic Galerkin technique formulated in section \ref{sec:SFEM}.
%
With PC characterized structural parameters, the probabilistic hysteretic behavior of restoring
force  versus  inter-story  drift  can  be  intrusively  modeled  following  the
stochastic Galerkin technique formulated in section \ref{sec:SFEM}.
%
Figure  \ref{figure_probabilisitc_constitutive_relation} shows the probabilistic
response of restoring force of the first floor under cyclic loading.
%
%
% \begin{figure}[!htbp]
% \centering
% \includegraphics[width=0.6\textwidth]{Figures/constitutive_relation_MC_verification.pdf}
% \caption{\label{figure_probabilisitc_constitutive_relation}  Intrusive inelastic
% probabilistic modeling of Armstrong-Frederick hysteretic behavior and  verification  with Monte Carlo simulation: Gaussian
% distributed  $Ha$  with  mean 1.76 $\times 10^{7} \ N/m$ and 15\% coefficient of
% variation (COV), $C_r = 17.6$.}
% \end{figure}
%
\begin{figure}[!htbp]
\centering
\subfloat[Uncertain $H_a$]{
\hspace{-0.8cm}
\includegraphics[width=0.53\textwidth]{Figures/constitutive_relation_uncertainHa_certainCr_MC_verification.pdf}}
\subfloat[Uncertain $H_a$ and $C_r$]{
\hspace{-0.2cm}
\includegraphics[width=0.53\textwidth]{Figures/constitutive_relation_uncertainHa_uncertainCr_MC_verification.pdf}}
\vspace{-2mm}
\caption{\label{figure_probabilisitc_constitutive_relation} Intrusive probabilistic modeling of Armstrong-Frederick hysteretic behavior and  verification  with Monte Carlo simulation: (a) Gaussian distributed  $Ha$  with  mean 1.76 $\times 10^{7} \ N/m$ and 15\% coefficient of variation (COV), $C_r = 17.6$. (b) Gaussian distributed  $Ha$  with  mean 1.76 $\times 10^{7} \ N/m$ and 15\% coefficient of variation (COV), Gaussian distributed $C_r$ with mean 17.6 and 15\% COV.}
\end{figure}
%
%
Verification  of  developed constitutive modeling is performed using 10,000 Monte
Carlo simulations and shown in Figure \ref{figure_probabilisitc_constitutive_relation}
as well.
%
It  can  be seen that PC-based intrusive probabilistic hysteresis modeling 
produces   almost  the  same  response  as  Monte  Carlo  simulations.
%
It  is  noted  that intrusive probabilistic approach is approximately 2000 times
faster than corresponding Monte Carlo modeling.
%
% From the quantiles plot shown in Figure \ref{figure_probabilisitc_constitutive_relation}(a),  


With  both  uncertain  seismic  motions (dimension 150) and uncertain structural
parameters  (dimension  4)  represented by Hermite PCs, probabilistic structural
displacement is described  in 154 dimensional probabilistic space of Hermite PCs.
%
The  unknown  time  varying  PC  coefficients, that contain all the information
about the probabilistic evolution of structural response, are intrusively solved
using developed Galerkin SEPFEM (section \ref{sec:SFEM}).
%
With  these  solved  PC  coefficients,  a polynomial chaos based surrogate model
is analytically established~\cite{Sudret2008}.
%
After  that, any probabilistic structural dynamic response can be
easily reconstructed.
%
Time evolving mean, standard deviation (SD) and correlation field
of any resulting field of interest
can   be   directly   evaluated   through  Equations  \ref{eq:PC_marginal_mean},
\ref{eq:PC_marginal_variance} and \ref{eq:PC_correlation_structure}.
%
By  efficiently  sampling  the  PC surrogate model, marginal or joint PDF of any
structural  response  of  interest  can  also be obtained through kernel density
estimation.



For  example,  Figure    \ref{figure_time_evolving_statistics}   shows  the time
evolving mean and standard deviation (SD) of the first and top floor deformation
relative to the ground.
%
%
\begin{figure}[!htbp]
\centering
\subfloat[Displacement mean]{
\includegraphics[width=0.55\textwidth]{Figures/mean_disp_longer.pdf}}
\subfloat[Displacement SD]{
\includegraphics[width=0.55\textwidth]{Figures/std_disp_longer.pdf}}
\caption{\label{figure_time_evolving_statistics} Time evolving mean and standard
deviation (SD) of the first and top floor deformation relative to the ground.}
\end{figure}
%
%
Due  to  inelastic, elastic-plastic response,  uncertain permanent deformation is observed
in both mean and standard deviation of floor deformation.
%
It is noted that the  deformation  of top floor presents much larger variability
than that of the first floor.




Two  typical  engineering demand parameters (EDPs) are selected for seismic risk
analysis:  Maximum  inter-story  drift  ratio (MIDR) and Peak floor acceleration
(PFA) \cite{Miranda2005,Miranda2006}.
%
We define MIDR as a function of probabilistic dynamic floor displacement:
%
\begin{equation}
MIDR_i(\theta) = \max_{t\in [0, T]} \Big\{\frac {|u_i(t, \theta) - u_{i-1}(t, \theta)|}{H_i} \Big\}
\label{eq:MIDR_floor}
\end{equation}
%
\begin{equation}
MIDR(\theta) = \max_{i\in [1, 8]}\max_{t\in [0, T]} \Big\{ \frac {|u_i(t, \theta) - u_{i-1}(t, \theta)|}{H_i} \Big\}
\end{equation}
%
%
%
where  $MIDR_i(\theta)$  and  $u_i(t,  \theta)$  are  the probabilistic MIDR and
displacement of the $i^{th}$ floor, respectively, and $H_i$ is the floor height,
while  probabilistic  MIDR  of  the  whole  shear  frame  structure  is given as
$MIDR(\theta)$.
%

Probabilistic floor accelerations are defined as :
%
%
%PFA is defined as a function of probabilistic dynamic floor acceleration: 
%
\begin{equation}
PFA_i(\theta) = \max_{t \in [0, T]} \{|\ddot{u}_i(t, \theta)|\}
\end{equation}
%
\begin{equation}
PFA(\theta) = \max_{i\in [1, 8]} \max_{t \in [0, T]} \{|\ddot{u}_i(t, \theta)|\}
\label{eq:PFA_overall}
\end{equation}
%
%
where  $PFA_i(\theta)$ and $\ddot{u}_i(t, \theta)$ are the probabilistic PFA and
acceleration  of  the  $i^{th}$  floor, respectively, while $PFA(\theta)$ is the
probabilistic PFA of the whole structure.
%
Since    both  probabilistic  displacements  $u_i(t,  \theta)$ and probabilistic
accelerations  $\ddot{u}_i(t,  \theta)$  are  well defined through resulting PC 
coefficients, probabilistic response of MIDR and PFA are readily
available through Equations \ref{eq:MIDR_floor} to \ref{eq:PFA_overall}.



%
For  example problem, the   probability  density  evolution  of  $MIDR(\theta)$ 
is shown in Figure \ref{figure_probability_evolution}.
%
\begin{figure}[!htbp]
\subfloat[PDF surface]{
\hspace{-1cm}
\includegraphics[width=0.65\textwidth]{Figures/MIDR_PDF_evolution.pdf}}
\subfloat[PDF contour]{
\includegraphics[width=0.50\textwidth]{Figures/Contour_PDF_evolution_MIDR.pdf}}
\caption{\label{figure_probability_evolution}  Time evolving probability density
function (PDF) of MIDR for frame structure.}
\end{figure}
%
%
At  $t=0s$,  the  structure is deterministically at rest, therefore,  the PDF of
$MIDR$ tends to infinity, i.e., a delta function centered at zero and as such is
not shown in Figure~\ref{figure_probability_evolution}.
% 
Figure  \ref{figure_evolving_PDF_combine}  shows typical PDFs at three different
times.
%
\begin{figure}[!htbp]
\centering
\includegraphics[width=0.6\textwidth]{Figures/PDF_MIDR_combine.pdf}
\caption{\label{figure_evolving_PDF_combine}  PDF  of  MIDR  at different times:
$t=15$s, $21$s and $29$s.}
\end{figure}
%
It  can  be  observed  that  PDF of MIDR is dispersing during  first half of the
seismic  loading,  while  toward the end of the loading, it shows high kurtosis,
due to reduced variation in input excitations.
%

The    PDFs    of  MIDR  of several  different  floors  (1st,  3rd,  6th  and top
floor)    and    the    whole    frame    structure  are   shown   in   Figure
\ref{figure_MIDR_different_stories}(a).
%
%
It  is  observed  that the mean of MIDR increases along with larger dispersion,
from the top to the bottom floor.
%
This  is  expected considering the increase of shear force from the top floor to
the base.
%
The  MIDR  PDF  of  the  first  floor  almost  overlaps  with  that of the whole
structure, which indicates that the maximum inter-story drift happens at first
floor.
%
From   the   probabilistic   distribution   of   MIDR,   exceeding   probability
$P(EDP>z|\boldsymbol \Gamma)$ can be obtained.
%
Combining  exceeding  probability  and  scenario  rate,  EDP  hazard  of MIDR is
calculated  using  Equation  \ref{eq:time_domain_risk_combine}  and  is shown in
Figure \ref{figure_MIDR_different_stories}(b).
%
%
\begin{figure}[!htbp]
\centering
\subfloat[PDF MIDR]{
\hspace*{-0.7cm}
\includegraphics[width=0.53\textwidth]{Figures/MIDR_distribution_different_floors.pdf}}
\subfloat[EDP hazard of MIDR]{
\hspace*{-0.7cm}
\includegraphics[width=0.53\textwidth]{Figures/Risk_MIDR.pdf}}
\caption{\label{figure_MIDR_different_stories} PDF and annual exceedance rate of
MIDR between different story over the whole loading history.}
\end{figure}
%
It can be seen that the demand of MIDR is dominantly controlled by lower floors,
e.g., the 1st and 3rd floor.
%



In  addition  to  PDFs of MIDR, PDFs  of  PFA for different floors and the whole
frame      structure      are      developed     and     shown     in     Figure
\ref{figure_PFA_different_stories}(a).
%
The  distributions of PFA of the 1st, 3rd and 6th floor are close to each other,
while the PFA of the top floor shows larger mean and variability.
%
The  PFA  distribution  of  the  top  floor  is  very close to that of the whole
structure,  which  indicates  the  top  floor  tends  to  experience the maximum
acceleration.
%
EDP hazard of PFA is shown in Figure \ref{figure_PFA_different_stories} (b).
%
The demand of PFA is dominantly controlled by the top floor.
%
\begin{figure}[!htbp]
\centering
\subfloat[PDF PFA]{
\hspace*{-0.7cm}
\includegraphics[width=0.53\textwidth]{Figures/PFA_distribution.pdf}}
\subfloat[EDP hazard of PFA]{
\hspace*{-0.7cm}
\includegraphics[width=0.53\textwidth]{Figures/Risk_PFA.pdf}}
\caption{\label{figure_PFA_different_stories}  PDF and annual exceedance rate of
PFA of different stories and the whole frame structure.}
\end{figure}





By assuming that  damage measure (DM) is a step function of EDP, seismic risk for damage
states  using  different  levels  of  MIDR  and  PFA  exceedance can be directly
determined from the EDP hazard curve.
%
As shown in Table \ref{tab_seismic_risk_summary}, seismic risk for MIDR$>1\%$ is
$3.83\times 10^{-3}$ and the risk for PFA$>1m/s^2$ is $1.92\times10^{-3}$.
%
\begin{table}[!htbp]
% \small
\footnotesize
% \scriptsize
\caption{\label{tab_seismic_risk_summary}Seismic   risk   of  damage  state  for
different levels of MIDR and PFA exceedance.}
\resizebox{0.98\hsize}{!}{
\begin{tabular}{cccccc}
%\hline
MIDR\textgreater{}0.5\% & MIDR\textgreater{}1\% & MIDR\textgreater{}2\% & PFA\textgreater{}0.5${\rm m/s^2}$ & PFA\textgreater{}1${\rm m/s^2}$ & PFA\textgreater{}1.5${\rm m/s^2}$ \\ 
\hline
6.66$\times 10^{-3}$                    & 3.83$\times 10^{-3}$                  & 9.97$\times 10^{-5}$                   & 6.65$\times 10^{-3}$                 & 1.92$\times 10^{-3}$               & 9.45$\times 10^{-5}$                 \\ 
%\hline
\end{tabular}}
\end{table}



As noted earlier, complete  probabilistic  structural  response,  including  both  marginal
distribution  and  correlation  information, is contained in PC coefficients, any
other  EDP  or other DM defined on multiple EDPs can also be developed with little
additional effort.
%
Figure  \ref{figure_joint_PFA_MIDR}  shows  the  2D  joint  PDF,  
$f(\text{MIDR}, \text{PFA}  |  \, \boldsymbol \Gamma)$ of two EDPs, 
MIDR and PFA, evaluated from
the PC-based surrogate model of probabilistic structure response.
%
%
\begin{figure}[!htbp]
\centering
\subfloat[Joint PDF of MIDR and PFA]{
\hspace*{-0.7cm}
\includegraphics[width=0.57\textwidth]{Figures/2D_EDP_PDF_1e5.pdf}}
\subfloat[Contour of 2D joint PDF]{
\hspace*{-0.7cm}
\includegraphics[width=0.53\textwidth]{Figures/2D_EDP_PDF_downview_1e5.pdf}}
\caption{\label{figure_joint_PFA_MIDR} 2D joint PDF of MIDR and PFA of the whole
shear frame structure.}
\end{figure}
%
It can be observed that in this case MIDR and PFA are positively correlated.
%
The correlation coefficient is $0.64$.
%





For   damage   measure   (DM)  defined  on  multiple  EDPs,  for  example,  $DM:
\{\text{MIDR}>z_1\, \vee\, \text{PFA}>z_2\}$, seismic risk can be evaluated as:

\begin{equation}
\lambda(\text{MIDR}>z_1 \vee \text{PFA}>z_2) 
= 
\overline{\lambda}\int_{\mathscr{D} } f(\text{MIDR}, \text{PFA} |\, \boldsymbol \Gamma) \, d \mathscr{D}  % \bigcup
\end{equation}
%
%
where  $\overline{\lambda}$  is  the  annual occurrence rate of seismic scenario
that   would   induce  ground  motion  population  $\boldsymbol  \Gamma$,  while
$\mathscr{D}$  is  the  integral  domain $(\text{MIDR}, \text{PFA}) \in $ $[z_1,
+\infty] \bigcup [z_2, +\infty]$ according to the definition of damage measure.
%
%
% Using such approach,   seismic   risk   for  damage  state  $DM:\{\text{MIDR}>1\%\,
% \vee\,\text{PFA}>1{\rm m/s^2}\}$  can be calculated  as  $4.20 \times 10^{-3}$, while the
% risk  for  damage state $DM:\{\text{MIDR}>1\%\, \wedge \,\text{PFA}>1{\rm m/s^2}\}$ is
% 60\% less, equal to $1.71 \times 10^{-3}$.
%

%  
%  {\huge \bf NOTE:
%  need to say in words what is 
%  $DM:\{\text{MIDR}>1\%\, \vee\,\text{PFA}>1{\rm m/s^2}\}$ is, and the other one too... 
% 
% so for example:
% 
% . . . 
% Damage state DM defined for MIDR that is greater than 1\% and  PFA that is greater than $1m/s^2$, 
% ($DM:\{\text{MIDR}>1\%\, \vee\,\text{PFA}>1{\rm m/s^2}\}$)
% . . . 
% 
% In a sense all mathematical equations and symbols have to be described in words if in text...
% 
% }

%%%% ===== Hexiang Reply: Has changed accordingly ==========
%
Using  such  approach,    seismic risk  for  damage  state DM defined for either
MIDR    greater    than    1\%    or    PFA   greater   than   $1{\rm m/s^2}$   (i.e.,
$DM:\{\text{MIDR}>1\%\,  \vee\,\text{PFA}>1{\rm m/s^2}\}$), can be calculated as
$4.20  \times  10^{-3}$,  while  the risk for damage state defined for both MIDR
greater  than  1\%  and PFA greater than $1{\rm m/s^2}$ (i.e., $DM:\{\text{MIDR}>1\%\,
\wedge  \,\text{PFA}>1{\rm  m/s^2}\}$) is  60\%  less,  equal  to  $1.71  \times
10^{-3}$.
%
Both of these risk values based on joint EDPs are rather different from the
ones calculated using single EDP.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}
\label{sec:conclusions}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A    time   domain  intrusive  probabilistic  seismic  risk  analysis  framework
for performance based earthquake engineering was described in some detail.
%
Methodology to simulate non-stationary stochastic seismic motions was 
presented.
%
The  presented  methodology is directly compatible with state-of-the-art seismic
source characterization.
%
Different  source,  path and site factors are explicitly accounted for in the
stochastic modeling of Fourier amplitude spectrum and Fourier phase derivative.
%
Both   uncertain   seismic  motions  and  uncertain  structural  parameters  are
characterized  as  random  process/field and represented with Hermite polynomial
chaos (PC) Karhunen-Lo{\`e}ve (KL) expansion.
% 
Direct  polynomial chaos based  Galerkin  intrusive  modeling  of 1D elastic-plastic response
was formulated  and  applied  to  simulate  the  uncertain hysteretic
behavior of restoring force versus inter-story drift for shear frame structure.
%
Formulations  for  random stiffness polynomial chaos coefficients were derived and incorporated
into stochastic Galerkin elastic-plastic finite element method.
%
% The  issue  of  PC characterization of random stiffness in current probabilistic
% elastic-plasticity based on Fokker-Planck-Kolmogorov (FPK) equation is resolved.
% %


Using  developed  stochastic elastic-plastic finite element method, probabilistic dynamic response
of  uncertain  structural  system  driven  by  uncertain  motions is intrusively
solved.
%
Following  that,  seismic  risk for damage measure defined on single or multiple
engineering demand parameter(s) was calculated.
%
The proposed framework is illustrated within seismic risk analysis of an eight-story
shear frame structure excited by uncertain strike-slip fault earthquakes.
%

Presented  new framework avoids the drawbacks of choosing and using intensity measure(s).
%
All  the  seismic  motion  characteristics and their uncertainties, for example,
uncertain peak ground acceleration (PGA), spectrum acceleration (Sa) and others, are captured by random process motions and
directly propagated into uncertain structural system.
%
Development   of ground motion prediction equations (GMPEs) for  potentially  new  intensity measures (IMs)  (e.g.,  Arias intensity or cumulative absolute velocity)  and repetitive  Monte Carlo   fragility simulations are
circumvented.
%
Though  most  of  current seismic risk analyses are performed for damage measure
defined  on  single  engineering demand parameter,  presented framework can also handle joint engineering demand parameters/failure
criteria without much additional effort.
%
It  is  found  that,  for  different  damage  measure  defined  on  joint  engineering demand parameters,
corresponding seismic risk significantly varies and is rather different from the
risk value for single engineering demand parameter.
%
Therefore,  considering  damage  measure  based  on  joint engineering demand parameters can be of great
interest in seismic risk analysis.
%
Future  work  will  focus  on  accuracy  and  efficiency  comparison between the
proposed framework and existing intensity measure based, non-intrusive seismic risk analysis and also applying the proposed framework to more realistic engineering structures. 
%


% Here's the list of references:
%
% \label{section:references}
\bibliography{refmech}
%

\end{document}
