PSHA with the OpenQuake engine
==============================

This chapter describes the mathematical framework for PSHA
implemented by the OpenQuake engine. More comprehensive descriptions
of Probabilistic Seismic Hazard Analysis can be found for example in
McGuire (2004) and USNRC (2012).

Basic concepts
--------------

Two main assumptions are at the base of all PSHA calculators included
in the OpenQuake engine: 

- seismicity in a region is described by a collection
  of *independent seismic sources* (i.e. the occurrence of an
  earthquake rupture in a source does not affect the probability of
  earthquake occurrence in the other sources)
- each source generates *independent earthquake ruptures* (i.e. the
  occurrence of an earth- quake rupture in a source does not affect the
  probability of occurrence of the other potential earthquake ruptures
  in the same source)

The Classical, Event-Based, and Disaggregation analysis requires the
definition of two main components: the *seismic source model*, that
is a collection of seismic sources describing the seismic activity in
a region of interest, and the *ground motion model*, that is a
mathematical relationship defining the probability distribution of a
ground motion parameter at a site given the occurrence of an
earthquake rupture.

The design of a seismic source model involves the specification of a
number of sources whose main parameters are the geometry,
constraining the earthquake rupture locations, and the
*magnitude-frequency distribution*, defining the average annual
occurrence rate over a magnitude range. A seismic source model
(*SSM*) can be therefore defined as a set of *I* seismic sources
(*Src*):

.. math::

 SSM = \{Src_1, Src_2, ..., Src_I\}

Chapter 3 provides a detailed description of the
different source typologies supported by the OpenQuake engine. However,
independently of the typology, in a PSHA each source undergoes a
discretization process which effectively generates a number of
distinct earthquake ruptures. A generic *i*-th source defines
therefore a set of *J* earthquake ruptures:

.. math::

 Src_i = \{Rup_{i1} ,Rup_{i2}, ..., Rup_{iJ}\}

Classical PSHA
--------------

The classical PSHA analysis allows calculating the probabilities of
exceeding, at least once in a given time span, and at a given site, a
set of ground motion parameter levels considering all possible
earthquake ruptures defined in a seismic source model. Such a list of
probability values is usually referred to as *hazard curve*.

We indicate with :math:`P(X \geq x|T)` the probability that a ground-motion
parameter *X* exceeds, at least once in a time span :math:`T` , a level
:math:`x`. :math:`P(X \geq x|T)` can be computed as 1 minus the probability that
none of the sources is causing a ground motion exceedance. By
assuming *independent sources*, the probability that none of the sources is 
causing an exceedance is equal to the product of the probabilities that each
source does not cause an exceedance, that is:

.. math::

 P(X \geq x|T) &= 1 - P_{src1}(X<x|T) \times P_{src2}(X<x|T) \times ... \times P_{srcI}(X<x|T)\\
               &= 1 - \prod^{I}_{i=1}{P_{src_i}(X<x|T)}

where :math:`P_{src_i}(X<x|T)` is the probability that the *i*-th
source is not causing an exceedance and :math:`I` is the total number of
sources in the source model.

By further assuming each source generates *independent earthquake
ruptures*, we can compute :math:`P_{src_i}(X<x|T)` as the product of
the probabilities that each rupture does not cause an exceedance,
that is:

.. math::

 P_{src_i}(X<x|T) &= P_{rup_{i1}}(X<x|T) \times P_{rup_{i2}}(X<x|T) \times ... \times P_{rup_{ij}}(X<x|T)\\
                  &= \prod^{J_i}_{j=1}{P_{rup_{ij}}(X<x|T)}

where :math:`P_{rup_{ij}}(X<x|T)` is the probability that the *j*-th
rupture in the *i*-th source is not causing an exceedance and :math:`J_i`
is the total number of ruptures generated by the *i*-th source.

Intuitively, the fact that a rupture does not cause any exceedance in
a given time span T can be due to the fact that the rupture does not
occur at all or that the rupture occurs once but without causing an
exceedance, or that the rupture occurs twice but both times without
causing an exceedance, and so on. Given that all these events are
mutually exclusive, by using the total probability theorem we can
write:

.. math::

 P_{rup_{ij}}(X<x|T) &= P_{rup_{ij}}(n=0|T)+P_{rup_{ij}}(n=1|T) \times P(X<x|rup_{ij}) + P_{rup_{ij}}(n=2|T) \times P(X<x|rup_{ij})^2 +...\\
                     &= \sum^{∞}_{k=0}{P_{rup_{ij}}(k|T) \times P(X<x|rup_{ij})^k}

where :math:`P_{rup_{ij}}(k|T)` is the probability that the *j*-th
rupture in the *i*-th source is occurring :math:`k` times in time span :math:`T`
and :math:`P(X<x|rup_{ij})` is the conditional probability that
parameter :math:`X` is not exceeding level :math:`x` given an occurrence of
:math:`rup_{ij}`.

Relying therefore on the assumptions of independent sources and
independent earthquake ruptures generated by each source, we can
compute the probability of at least one ground motion exceedance as:

.. math::

 P(X \geq x|T) &= 1-\prod^{I}_{i=1}\prod^{J_i}_{j=1}{P_{rup_{ij}}(X<x|T)}\\
               &= 1-\prod^{I}_{i=1}\prod^{J_i}_{j=1}\sum^{∞}_{k=0}{P_{rup_{ij}}(k|T) \times P(X<x|rup_{ij})^k}

It is worth noticing how equation 2.6 allows the
calculation of hazard curves for any temporal occurrence model. More
precisely, only the probability of the number of rupture occurrences
:math:`P_{rup_{ij}}(k|T)` depends on the temporal occurrence model, while
the procedure for calculating the final probability of ground motion
exceedance remains the same. This also means that is possible to
define in the same source model sources with different temporal
occurrence models (for instance time-dependent and time-independent
sources).

Poissonian source model
***********************

If we now assume a source model to consist only of Poissonian
sources, we can write :math:`P_{rup_{ij}}(k|T)` for every rupture as:

.. math::

 P_{rup_{ij}}(k|T) = e^{-\nu_{ij}T}\frac{(\nu_{ij}T)^k}{k!}

where :math:`nu_{ij}` is the average annual occurrence rate for the *j*-th
rupture in the *i*-th source. We can then place equation
above and thus write:

.. math::

 P_{rup_{ij}}(X<x|T) &= \sum^{∞}_{k=0}{e^{-\nu_{ij}T}\frac{(\nu_{ij}T)^k}{k!}} \times P(X<x|rup_{ij})^k\\
                     &= e^{-\nu_{ij}T} \sum^{∞}_{k=0}{\frac{(\nu_{ij}T \times P(X<x|rup_{ij}))^k}{k!}}

Making use of the property:

.. math::

 e^x = \sum^{∞}_{k=0}{\frac{x^k}{k!}}

we can rewrite 2.8 as:

.. math::

 P_{rup_{ij}}(X<x|T) &= e^{-\nu_{ij}T}e^{\nu_{ij}T \times P(X<x|rup_{ij})}\\
                     &= e^{-\nu_{ij}T \times (1-P(X<x|rup_{ij}))}\\
                     &= e^{-\nu_{ij}T \times P(X \geq x|rup_{ij})}

By now recognizing that, according to the Poissonian distribution,
the probability of at least one occurrence (that is one or more) in a
time span :math:`T` of :math:`rup_{ij}` is:

.. math::

 P_{rup_{ij}}(n \geq 1|T) = 1-e^{-\nu_{ij}T}

we can write equation 2.10 as:

.. math::

 P_{rup_{ij}}(X<x|T) = (1-P_{rup_{ij}}(n \geq 1|T))^{P(X \geq x|rup_{ij})}

By placing equation 2.12 in 2.6, we can write:

.. math::

 P(X \geq x|T) = 1-\prod^{I}_{i=1}\prod^{J_i}_{j=1}(1-P_{rup_{ij}}(n \geq 1|T))^{P(X \geq x|rup_{ij})}

Equation 2.13 is used by the OpenQuake engine for the
calculation of hazard curves when performing Classical PSHA with a
Poissonian source model. To our knowledge, this equation has been
first proposed by Field et al. (2003), derived
from the traditional rate-based formulation converted in terms of
probabilities (their equation A8). Instead, we derive it from the
assumptions of a source model consisting of independent sources,
independent earthquake ruptures generated by each source, and
ruptures obeying to a Poissonian temporal occurrence model.

Equivalence with the rate-based equation
****************************************

It is worth noticing how equation 2.13 is
equivalent to the more traditional rate-based hazard equation
(McGuire, 1995). Indeed, by assuming ground
motion occurrence to follow a Poissonian distribution in time, and
indicating with *n* the mean annual rate of exceeding a ground motion
level x, we can write:

.. math::

 P(X \geq x|T) = 1-e^{-\nu T}

We can also rewrite 2.13 as:

.. math::

 P(X \geq x|T) &= 1-\prod^{I}_{i=1}\prod^{J_i}_{j=1}(1-P_{rup_{ij}}(n \geq 1|T))^{P(X \geq x|rup_{ij})}\\
               &= 1-\prod^{I}_{i=1}\prod^{J_i}_{j=1}e^{-\nu_{ij}T \times P(X \geq x|rup_{ij})}\\
               &= 1-e^{-\sum^{I}_{i=1}\sum^{J_i}_{i=1}\nu_{ij}T \times P(X \geq x|rup_{ij})}

The equivalence between equations and 2.15 is possible if and only if:

.. math::

 \nu = \sum^{I}_{i=1}\sum^{J_i}_{i=1}\nu_{ij} \times P(X \geq x|rup_{ij})

Assuming now, for the sake of simplicity, that a rupture is
completely characterized by magnitude and distance from a site, we
can write the rate of occurrence of the *j*-th rupture as:

.. math::

 \nu_{ij} = \nu_{i} \times f_{i}(m,r)

where :math:`nu_{i}` is the total occurrence rate for the *i*-th source, and
:math:`f_i(m,r)` is the probability, for the *i*-th source, of
generating a rupture of magnitude :math:`m` and distance :math:`r`. By placing
2.17 in 2.16 and by replacing
the discrete summation over ruptures with a continuous integral over
magnitude and distance, we can write:

.. math::

 \nu = \sum^{I}_{i=1} \nu_i \int\int f_i(m,r)P(X \geq x|m,r)dmdr

which is the traditional equation for calculating ground motion
exceedance rates (McGuire, 1995).

Event-based PSHA
----------------

The goal of an Event-based PSHA is to simulate seismicity in a region
as described by a source model and to simulate ground shaking on a
set of locations accordingly with a ground motion model. In both
cases, simulation involves a Monte Carlo (i.e. random) sampling
procedure.

Seismicity is simulated by generating a *stochastic event set* (also
known as *synthetic catalog*) for a given time span :math:`T`. For each
rupture generated by a source, the number of occurrences in a time
span :math:`T` is simulated by sampling the corresponding probability
distribution as given by :math:`P_{rup}(k|T)`. A stochastic event set is 
therefore a *sample* of the full population of ruptures as defined by a 
seismic source model. Each rupture is present zero, one or more times, 
depending on its probability. Symbolically, we can define a stochastic event
set :math:`SES` as:

.. math::

 SES(T) = \{k \times rup,\ k\ \sim\ P_{rup}(k|T)\ \forall\ rup\ in\ Src\ \forall\ Src\ in\ SSM\}

where :math:`k`, the number of occurrences, is a random sample of
:math:`P_{rup}(k|T)`, and :math:`k \times rup` means that rupture 
:math:`rup` is repeated :math:`k` times in the stochastic event set.

Given an earthquake rupture, the simulation of ground shaking values
on a set of locations :math:`x=(x_1,x_2,.....,x_N)` 
forms a *ground motion field*. In a Event-based PSHA, for each rupture in a
stochastic event set, the ground motion field is obtained by sampling
the probability distribution defined by the ground motion model. As
described in Chapter 4, the ground motion
distribution at a site is assumed to be a Normal distribution. The
aleatory variability is described in terms of an *inter-event* (also
known as *between-events*) standard deviation (:math:`\tau`) and intra-event
(also known as *within-event*) standard deviation (:math:`\sigma`). The
simulation of a ground motion field is therefore the result of the
summation of three terms, the logarithmic mean of the ground motion
distribution:

.. math::

 \mu = (\mu_1,\mu_2,...,\mu_N)

the inter-event variability:

.. math::

 \eta = (\eta,\eta,...,\eta),\ where\ \eta \sim N(0,\tau)

and the intra-event variability:

.. math::

 \epsilon = (\epsilon_1,\epsilon_2,...,\epsilon_N)\ \sim N(0,\Sigma)
   
where:

.. math::

 \Sigma = \left [ \begin{matrix} \sigma_1^2 & \sigma_1\sigma_2\rho_{12} & ... & \sigma_1\sigma_N\rho_{1N} \\ 
                        \sigma_2\sigma_1\rho_{21} & \sigma_2^2 & ... & \sigma_2\sigma_N\rho_{2N} \\
                        ... & ... & ... & ... \\
                        \sigma_N\sigma_1\rho_{N1} & \sigma_N\sigma_2\rho_{N2} & ... & \sigma_N^2 \end{matrix}\right]

It is worth noticing how the inter-event variability, uniform for all
sites given an earthquake rupture, is drawn from an univariate normal
distribution of mean 0 and standard deviation *t*, while the
intra-event variability is a random sample of a multivariate normal
distribution of mean :math:`0` and covariance matrix :math:`\Sigma`. 
:math:`\Sigma` is a
diagonal matrix when considering no correlation in the intra-event
variability, while it has non-zero off-diagonal elements when
considering a correlation model (:math:`\rho`) for the intra-event aleatory
variability.

Calculation of hazard curves from ground motion fields
******************************************************

The ground motion fields simulated for each rupture in a stochastic
event set can be used to compute hazard curves. Indeed, indicating
with :math:`T_0` the duration associated with a stochastic event
set, and with :math:`K` the number of ground motion fields (and associated
ruptures) simulated in time :math:`T_0`, we can compute the rate of
exceedance of a ground motion level :math:`x` at a site as (Ebel and Kafka,
1999):

.. math::

 \nu = \frac{\sum^{K}_{k=1}H(x_k-x)}{T_0}

where :math:`H` is the Heaviside function and :math:`x_k` is the ground motion
parameter value at the considered site associated with the *k*-th
ground motion field. The exceedance rate obtained from equation 2.24 can then be 
used to compute the probability of at least one occurrence in any time span 
:math:`T` , accordingly with the Poissonian distribution, using equation 2.14.
This approach is possible when the source model from which the
stochastic event set is generated is Poissonian.

As the stochastic event set duration :math:`T_0` increases,
equation 2.24 provides an increasingly more
accurate estimate of the actual rate of exceedance. A larger
:math:`T_0` can be achieved not only by simulating a single
stochastic event set with longer duration, but also by simulating
multiple stochastic event sets. These can then be joined together to
form a stochastic event which has a large enough duration to provide
a stable estimate of the rates of ground motion exceedance.

Disaggregation
--------------

The disaggregation analysis allows investigating how the different
earthquake ruptures defined in a source model contribute to the
probability of exceeding a certain ground motion level :math:`x` at a given
site. Given the very large number of earthquake ruptures associated
with a source model, contributions cannot be investigated on a
rupture by rupture basis but a classification scheme is used instead.
Ruptures are classified in terms of the following parameters:

- magnitude (:math:`M`)
- distance to rupture surface-projection (Joyner-Boore distance) (:math:`r_{jb}`)
- longitude and latitude of rupture surface-projection closest
  point (:math:`\lambda` , :math:`\phi`) 
- tectonic region type (:math:`TRT`)

For each earthquake rupture, the associated conditional probability
of ground motion exceedance – :math:`P(X \geq x|T,rup)` – is computed for different 
:math:`\epsilon` bins, where :math:`\epsilon` is the difference, in terms of number of total standard 
deviations, between :math:`x` and the mean ground motion :math:`\mu` as predicted by the ground
motion model, that is:

.. math::

 \epsilon = \frac{x-\mu}{\sigma_{total}}

The disaggregation in terms of *e* allows investigating how the
different regions of the ground motion distributions contribute to
the probability of exceedance.

The rupture parameters (:math:`M`, :math:`r_{jb}`, :math:`\lambda`, :math:`\phi`
, :math:`TRT`) together with the :math:`\epsilon` parameter effectively create
a 6-dimensional model space which, discretized into a number of bins,
is used to classify the probability of exceedance for different
combination of rupture parameters and :math:`\epsilon` values.

For a given model space bin :math:`m` = (:math:`M`, :math:`r_{jb}`, :math:`\lambda`, :math:`\phi`
, :math:`TRT`) the probability of exceeding level :math:`x`
at least once in a time span :math:`T` is computed using equation
2.6, that is:

.. math::

 P(X>x|T,m) = 1- \prod^{I}_{i=1}\prod^{J_i}_{j=1} 
  \begin{cases} P_{rup_{ij}}(X<x|T)\ & \text {if } rup_{ij}\ \in\ m \\ 
                                  1\ & \text {if } rup_{ij}\ \notin\ m \end{cases}

In other words, if a rupture belongs to the considered bin, then the
probability of not causing a ground motion exceedance is computed
according to equation 2.6, otherwise the
probability is 1 (that is, given that the rupture does not belong to
the bin it can never cause a ground motion exceedance).

Disaggregation histograms
*************************

Disaggregation values as given by equation 2.26
can be aggregated in order to investigate earthquake rupture
contributions over a reduced model space. The following
disaggregation histograms are provided by the OpenQuake engine.

Magnitude disaggregation:

.. math::

 P(X>x|T,M) = 1 - \prod_{r_{jb}}\prod_{\lambda}\prod_{\phi}\prod_{TRT}\prod_{\epsilon}(1-P(X>x|T,m))

Distance disaggregation:

.. math::

 P(X>x|T,r_{jb}) = 1 - \prod_{M}\prod_{\lambda}\prod_{\phi}\prod_{TRT}\prod_{\epsilon}(1-P(X>x|T,m))

Tectonic region type disaggregation:

.. math::

 P(X>x|T,TRT) = 1 - \prod_{M}\prod_{r_{jb}}\prod_{\lambda}\prod_{\phi}\prod_{\epsilon}(1-P(X>x|T,m))

Magnitude-Distance disaggregation:

.. math::

 P(X>x|T,M,r_{jb}) = 1 - \prod_{\lambda}\prod_{\phi}\prod_{TRT}\prod_{\epsilon}(1-P(X>x|T,m))

Magnitude-Distance-Epsilon disaggregation:

.. math::

 P(X>x|T,M,r_{jb},\epsilon) = 1 - \prod_{\lambda}\prod_{\phi}\prod_{TRT}(1-P(X>x|T,m))

Longitude-Latitude disaggregation:

.. math::

 P(X>x|T,\lambda,\phi) = 1 - \prod_{M}\prod_{r_{jb}}\prod_{TRT}\prod_{\epsilon}(1-P(X>x|T,m))

Longitude-Latitude-Magnitude disaggregation:

.. math::

 P(X>x|T,\lambda,\phi,M) = 1 - \prod_{r_{jb}}\prod_{TRT}\prod_{\epsilon}(1-P(X>x|T,m))

Longitude-Latitude-Tectonic Region Type disaggregation:

.. math::

 P(X>x|T,\lambda,\phi,TRT) = 1 - \prod_{M}\prod_{r_{jb}}\prod_{\epsilon}(1-P(X>x|T,m))

All the above equations are based on the assumption that earthquake
ruptures in different bins are independent, therefore probabilities
can be aggregated by using the multiplication rule for independent
events. The probability of a ground motion exceedance over a reduced
model space is computed simply as 1 minus the probabilty of
non-exceedance over the remaining model space dimensions.

Comparison between OpenQuake engine disaggregation and *traditional* disaggregation
***********************************************************************************

The traditional disaggregation analysis as commonly known in
literature (e.g. Bazzurro and Cornell, 1999)
differs from the one provided by the OpenQuake engine. Indeed, a
disaggregation analysis tipically provides the conditional
probability of observing an earthquake scenario of given properties
(magnitude, distance, epsilon, ...) given that a ground motion
exceedance is occurred, which can be written (following the same
notation used in this chapter) as:

.. math::

 P(m|X>x)

On the contrary, the OpenQuake engine (as described in equation
2.26) provides the conditional probability of
observing at least one ground motion exceedance in a time span :math:`T`
given the occurrence of earthquake ruptures of given properties
:math:`m`, that is:

.. math::

 P(X>x|T,m)

The probabilities given in equations 2.35 and
2.36 are clearly different. Indeed, for different
:math:`m`, values given by equation 2.35 must sum up
to 1, while this is not the case for equation
2.36. For the former equation different :math:`m`
represent mutually exclusive events, while for the latter they
represent independent events.

When considering a Poissonian source model it is possible however to
derive equation 2.35 from equation 2.36. Indeed,
indicating with :math:`\nu_m` the rate of ground motion exceedance
:math:`(X>x)` associated with earthquake ruptures of properties :math:`m` and
with :math:`\nu` the rate of ground motion exceedance associated with all
earthquake ruptures, we can write equation 2.35 as:

.. math::

 P(m|X>x)=\frac{\nu_m}{\nu}

By solving the Poissonian equation 2.14 for the
rate of exceedance, we can write :math:`\nu_m` as:

.. math::

 \nu_m=-\frac{\ln(1-P(X>x|T,m))}{T}

:math:`\nu` can be obtained using the same equation above but considering
:math:`P(X>x|T)` instead of :math:`P(X>x|T,m)`, where :math:`P(X>x|T)` 
is obtained by aggregating, using the multiplicative rule,
the probabilities over the different :math:`m`, that is:

.. math::

 P(X>x|T)=1-\prod_{m}(1-P(X>x|T,m))

By computing :math:`\nu_m` and :math:`\nu` from :math:`P(X>x|T,m)` it is
hence possible to obtain the more traditional disaggregation results
as given in equation 2.35.