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):
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:
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
where
By further assuming each source generates independent earthquake
ruptures, we can compute
where
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:
where
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:
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
Poissonian source model#
If we now assume a source model to consist only of Poissonian
sources, we can write
where
Making use of the property:
we can rewrite 2.8 as:
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
we can write equation 2.10 as:
By placing equation 2.12 in 2.6, we can write:
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:
We can also rewrite 2.13 as:
The equivalence between equations and 2.15 is possible if and only if:
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:
where
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
where
Given an earthquake rupture, the simulation of ground shaking values
on a set of locations
the inter-event variability:
and the intra-event variability:
where:
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
Rupture sampling: how does it work?#
In this section we explain how the sampling of ruptures in event based calculations works, at least for the case of Poissonian sources. As an example, consider the following point source:
>>> from openquake.hazardlib import nrml
>>> src = nrml.get('''\
... <pointSource id="1" name="Point Source"
... tectonicRegion="Active Shallow Crust">
... <pointGeometry>
... <gml:Point><gml:pos>179.5 0</gml:pos></gml:Point>
... <upperSeismoDepth>0</upperSeismoDepth>
... <lowerSeismoDepth>10</lowerSeismoDepth>
... </pointGeometry>
... <magScaleRel>WC1994</magScaleRel>
... <ruptAspectRatio>1.5</ruptAspectRatio>
... <truncGutenbergRichterMFD aValue="3" bValue="1" minMag="5" maxMag="7"/>
... <nodalPlaneDist>
... <nodalPlane dip="30" probability="1" strike="45" rake="90" />
... </nodalPlaneDist>
... <hypoDepthDist>
... <hypoDepth depth="4" probability="1"/>
... </hypoDepthDist>
...</pointSource>''', investigation_time=1, width_of_mfd_bin=1.0)
The source here is particularly simple, with only one seismogenic depth and one nodal plane. It generates two ruptures,
because with a width_of_mfd_bin
of 1 there are only two magnitudes in the range from 5 to 7:
>>> [(mag1, rate1), (mag2, rate2)] = src.get_annual_occurrence_rates()
>>> mag1
5.5
>>> mag2
6.5
The occurrence rates are respectively 0.009 and 0.0009. So, if we set the number of stochastic event sets to 1,000,000:
>>> num_ses = 1_000_000
we would expect the first rupture (the one with magnitude 5.5) to occur around 9,000 times and the second rupture (the one with magnitude 6.5) to occur around 900 times. Clearly the exact numbers will depend on the stochastic seed; if we set:
>>> np.random.seed(42)
then we will have (for investigation_time = 1
):
>>> np.random.poisson(rate1 * num_ses * 1)
8966
>>> np.random.poisson(rate2 * num_ses * 1)
921
These are the number of occurrences of each rupture in the effective investigation time, i.e. the investigation time multiplied by the number of stochastic event sets and the number of realizations (here we assumed 1 realization).
The total number of events generated by the source will be number_of_events = sum(n_occ for each rupture)
i.e. 8,966 + 921 = 9,887, with ~91% of the events associated to the first rupture and ~9% of the events associated to
the second rupture.
Since the details of the seed algorithm can change with updates to the the engine, if you run an event based calculation
with the same parameters with different versions of the engine, you may not get exactly the same number of events, but
something close given a reasonably long effective investigation time. After running the calculation, inside the
datastore, in the ruptures
dataset you will find the two ruptures, their occurrence rates and their integer number of
occurrences (n_occ
). If the effective investigation time is large enough then the relation n_occ ~ occurrence_rate * eff_investigation_time
will hold. If the effective investigation time is not large enough, or the occurrence rate is extremely small, then you
should expect to see larger differences between the expected number of occurrences and n_occ
, as well as a strong
seed dependency.
It is important to notice than in order to determine the effective investigation time, the engine takes into account also
the ground motion logic tree and the correct formula to use is eff_investigation_time = investigation_time * num_ses * num_rlzs
where num_rlzs
is the number of realizations in the ground motion logic tree.
Just to be concrete, if you run a calculation with the same parameters as described before, but with two GMPEs instead
of one (and number_of_logic_tree_samples = 0
), then the total number of paths admitted by the logic tree will be 2
and you should expect to get about twice the number of occurrences for each rupture. Users wanting to know the
nitty-gritty details should look at the code, inside hazardlib/source/base.py, in the method
src.sample_ruptures(eff_num_ses, ses_seed)
.
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
where
As the stochastic event set duration
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
magnitude (
)distance to rupture surface-projection (Joyner-Boore distance) (
)longitude and latitude of rupture surface-projection closest point (
, )tectonic region type (
)
For each earthquake rupture, the associated conditional probability
of ground motion exceedance –
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 (
For a given model space bin
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:
Distance disaggregation:
Tectonic region type disaggregation:
Magnitude-Distance disaggregation:
Magnitude-Distance-Epsilon disaggregation:
Longitude-Latitude disaggregation:
Longitude-Latitude-Magnitude disaggregation:
Longitude-Latitude-Tectonic Region Type disaggregation:
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:
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
The probabilities given in equations 2.35 and
2.36 are clearly different. Indeed, for different
When considering a Poissonian source model it is possible however to
derive equation 2.35 from equation 2.36. Indeed,
indicating with
By solving the Poissonian equation 2.14 for the
rate of exceedance, we can write
By computing