Theme V – Models and Techniques forAnalyzing Seismicity
Basic models of seismicity: Spatiotemporal models

Jiancang Zhuang1 • Maximilian J. Werner2 • Sebastian Hainzl3 •David Harte4 • Shiyoung Zhou5

1.
2.
3.
4.
5.

Institute of Statistical Mathematics
Department of Geosciences, Princeton University
Institute of Geosciences, University of Potsdam
Statistics Research Associates
Department of Geophysics, Peking University

How to cite this article:
Zhuang, J., M.J. Werner, S. Hainzl, D. Harte, and S. Zhou (2011), Basic models ofseismicity: spatiotemporal models, Community Online Resource for Statistical
Seismicity Analysis, doi:10.5078/corssa-07487583. Available at
http://www.corssa.org.
Document Information:
Issue date: 26 July 2011 Version: 1.0

2

www.corssa.org

Contents
1
2
3
4
5
6
7

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Formulation of the spatiotemporal ETAS models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Technical points related to the spatiotemporal ETAS model . . . . . . . . . . . . . . . . . . . . . . . . . .
Stochastic reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Further reading on the ETAS model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The EEPAS model and the double branching model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3
5
7
11
14
16
17

Basic models of seismicity: Spatiotemporal models

3

Abstract In this article, we present a review of spatiotemporal point-process models, including the epidemic type aftershock sequence (ETAS) model, the EEPAS(Every Earthquake is Precursor According to Scale) model, the double branchingmodel, and related techniques. Here we emphasize the ETAS model, because it has
been well studied and is currently a standard model for testing hypotheses related
to seismic activity.
1 Introduction
We assume that the reader has read the articles in this series on purely temporalmodels and the time independent spatial models (nonhomogeneous Poisson models) and that the reader is familiar with related concepts of point process models,
especially the concepts of conditional intensity, likelihood, temporal ETAS model,
and methods for estimating spatial rates of nonhomogeneous Poisson models. In thisarticle, we continue with spatiotemporal models.
There have not been very many spatiotemporal models to describe seismicity.
One reason is that heavy and complicated numerical computations are always involved in the implementation. Among the few spatiotemporal models that havebeen proposed, only the spatiotemporal epidemic-type aftershock sequence (ETAS)
model has been extensively studied in the context of earthquake short-term clustering. Mathematically, the model is in the class of marked branching processes withimmigration. Even though its formulation is based on empirical statistics, recent
studies by Console et al. (2007) and Iwata (2010) suggest that the ETAS model isthe best model for describing short-term seismicity, outperforming stochastic implementations of the physics-based model of Dieterich (1994). Both the EEPAS(“Every Earthquake is a Precursor According to Scale”, Rhoades and Evison 2004)
and the double branching models (Marzocchi and Lombardi 2008) have aspects incommon with the ETAS model, i.e., using the ETAS model as the standard clustering component. For these reasons, we focus on issues related to the spatiotemporal
ETAS model, such as the estimation of the background rate, stochastic declustering,stochastic reconstruction, etc.
The temporal ETAS model was suggested by Ogata (1988) with the conditional
intensity function of
λ(t) = µ + K0

X exp[α(m − m0 )]ti <t(t − ti + c)p,(1)

where µ (shocks per unit time) represents the rate of background seismicity, thesummation is taken over aftershocks occurring before time t, and m0 represents
the cut-off magnitude of the fitted data. In the above equation, the coefficient α

4www.corssa.org(magnitude−1 ) is a measure of the efficiency of a shock in generating aftershock
activity relative to its magnitude, K0 represents the productivity of an event ofthreshold magnitude m0 , and c (unit of time) and p are the parameters in the
Omori-Utsu law for describing the decay of the aftershock sequence.
Before Ogata generalized the temporal ETAS model to a space-time version,
Musmeci and Vere-Jones (1992) used space-time diffusion clustering models to analyze seismicity in Italy. The conditional intensity functions for these models havethe common formXλ(t, x, y) = µ(x, y) +gφ (t − ti , x − xi , y − yi , mi ),(2)
i:ti <t

where



1 x2
eαmi e−c t
y2
exp −
,
gφ (t, x, y, m) = A+
2πσx σy t2t σx2 σy2(3)

Aeαm e−ct t2 Cx Cy
gφ (t, x, y, m) = 2 2,
π (x + t2 Cx2 )(y 2 + t2 Cy2 )(4)

and A, α, σx , σy , Cx , and Cy are constants. For a fixed point (x, y), when t → ∞ theaftershocks in (3) and (4) decay with time according to t−1 e−ct and t−2 e−ct . Kagan(1991) and Rathbun (1993) also discussed slightly different forms.
The space-time ETAS model defined by Ogata (1998), now the generally accepted
definition, had the same general form as in (2), but where gφ (t, x, y, m) was defineddifferently, as
gφ (t, x, y, m) = κ(m)g(t)f (x, y|m).
(5)
In the above equation,
κ(m) = Aeα(m−m0 )(6)

is the expected number of aftershocks generated from a mainshock of magnitudem,p−1
tg(t) =(1 + )−p(7)cc
is probability density function of the lagged time distribution of aftershocks, and 2

1
x + y2f (x, y|m) =
f
(8)
πσ(m)
σ(m)
is the density function of the aftershock locations from a mainshock at the origin
with magnitude m.

Basic models of seismicity: Spatiotemporal models

5

2 Formulation of the spatiotemporal ETAS models
Even though there are some slight differences between the forms adopted by different
researchers, the common forms of these models can be outlined as follows:(a) The background events are regarded as the “immigrants” in the branching process of earthquake occurrence; their occurrence rate is assumed to be a functionof spatial location and magnitude, but not of time.
(b) Each event produces offspring events independently of the others. The expected
number of direct offspring from an individual event is assumed to depend on its
magnitude m, and will be denoted by κ(m), which is called its productivity.
(c) The probability distribution of the time until the appearance of a child event is afunction of the time lag from its parent, and is independent of magnitude; thus, its
probability density function is assumed to have the form g(t|τ ) = g(t − τ ) whereτ is the occurrence time of the ancestor. Moreover, the function is independentof what happens between τ and t.
(d) The probability distributions of the location (x, y) and magnitude m of a child
event are dependent on the magnitude m∗ and the location (ξ, η) of its parent.
These probability density functions are denoted by f (x−ξ, y−η|m∗ ) and s(m|m∗ ),respectively, where ξ, η and m∗ are the location and magnitude of the ancestor.
(e) The magnitudes of all the events, including background events and their offspring, are independent random variables drawn from the same probability distribution of density s(m).
In general, this class of marked branching point processes for earthquake occurrences can be represented by its conditional intensity function, which is written
Pr {an event in [t, t + dt) × [x, x + dx) × [y, y + dy) × [m, m + dm)| Ht }= λ(t, x, y, m) dt dx dy dm + o(dt dx dy dm),
(9)
where Ht denotes the space-time-magnitude occurrence history of the earthquakes
up to time t.
Based on the assumptions (a)–(e), the conditional intensity function for the spacetime model can be written asX
λ(t, x, y, m) = s(m) µ(x, y) +κ(mk )g(t − tk )f (x − xk , y − yk ; mk ) . (10){k: tk <t}In the above equation, µ(x, y) is the background intensity function, which is assumed to be independent of time. The functions g(t), f (x, y; mk ) and s(m) arerespectively the normalized response functions (i.e., probability density functions)of the occurrence time, the location, the magnitude of an offspring from an ancestor

6

www.corssa.org

of magnitude mk . And from the fact that the k-th event excites a non-stationary
Poisson process with intensity function κ(mk )g(t − tk )f (x − xk , y − yk ; mk )s(m), wenote that κ(mk ) is the expected number of children from a parent of size mk . Thisimplies that the number of children is a Poisson random variable with a mean ofκ(mk ). Note also that the probability function g(t) is independent of the magnitudeof the parent, as mentioned in assumption (c). Another explanation of (10) is thatthe risk of earthquake occurrence at time t and location (x, y) consists of the contribution from the background rate µ and the contributions from each previous event,ξ(t, x, y; ti , xi , yi ).
According to (10), the conditional intensity function for the model can be decomposed asλ(t, x, y, m) = s(m)λ(t, x, y),
(11)
where
s(m) = β e−β(m−m0 )(12)represents the Gutenberg-Richter law for earthquake magnitudes of m0 or larger inthe form of a probability density function, β is linked with the so-called b-value by
β = b ln 10, andX
λ(t, x, y) = µ(x, y) +κ(mk )g(t − tk )f (x − xk , y − yk ; mk ).
(13){k: tk <t}

In applications, the following specific functions are often used. The productivity
law is given by (6) and the Omori-Utsu law by (7).
For the spatial component, f (x, y; m), the following different functions have beenused:Model 1 [Rathbun 1993; Console et al. 2003]
1 − x2 +y2 2
e 2D ;
2πD2Model 2 [Ogata 1998; Zhuang et al. 2002]
f (x, y; m) =2 +y 21− 2 xα(m−m
0) ,
2D e
e
2πD2 eα(m−m0 )where the parameter α is the same one as in (6);Model 3 [Ogata 1998; Console et al. 2003]

−q
q−1
x2 + y 2f (x, y; m) =
1+
;
πD2
D2

f (x, y; m) =Model 4 [Ogata 1998; Zhuang et al. 2002, 2004]−qx2 + y 2
q−1
1 + 2 α(m−m0 );f (x, y; m) =πD2 eα(m−m0 )
D e

(14)

(15)

(16)

(17)

Basic models of seismicity: Spatiotemporal models

andModel 5 [Zhuang et al. 2005; Ogata and Zhuang 2006]−qx2 + y 2
q−1
1 + 2 γ(m−m0 )f (x, y; m) =.
πD2 eγ(m−m0 )
D e

7

(18)

The difference between Models 4 and 5 is that α is identical to the parameter α
in the productivity function κ(m), whereas γ is completely independent of κ(m).
As shown in Zhuang et al. (2004); Zhuang (2006), and Ogata and Zhuang (2006),
Model 5 (18) usually fits earthquake data the best among those five forms. However,
the other four forms are also frequently used by many researchers.
3 Technical points related to the spatiotemporal ETAS model
3.1 Maximum likelihood procedure
Given the background rate µ(x, y) and the observed earthquake catalog, the maximum likelihood estimates (MLE) of the model, θ̂ = (ν̂, Â, α̂, ĉ, p̂, D̂), are calculated
by maximizing the log-likelihood function
Z T ZZX
log L(θ) =log λθ (ti , xi , yi ) −λθ (t, x, y) dx dy dt,(19)
i:ti ∈[0,T ],(xi ,yi )∈S

0Swhere the index i runs over all the events occurring in the study region S andthe study time interval [0, T ]. The computational details can be found in Ogata(1998). Veen and Schoenberg (2008) developed an EM (expectation-maximization)
algorithm for a fast estimation by maximizing the expected log-likelihood.
In (19), the space-time window [0, T ] × S is called the target window, and theevents inside are called target events. Different from the target events, the index
k in (10) runs over all the recorded events in the catalog, whose range is called theauxiliary or complementary window. Usually the complementary window should betaken as large as possible. A more detailed discussion of the influence of the sizeof the complementary window on the estimated model parameters was provided by
Wang et al. (2010).
3.2 Thinning procedure
An interesting question in seismicity analysis is to what extent a given earthquake is
triggered by a previous earthquake. The so-called thinning procedure (see, e.g., Lewis
and Shedler 1979; Ogata 1981; Daley and Vere-Jones 2003), is an important tool forsuch an analysis and will be briefly described here. Below, we use the concept for8

www.corssa.org

assigning a “degree of parenthood” to earthquake pairs and to estimate the spatial
background density in the ETAS model. The proportion of the contribution fromevent i at the occurrence of (tj , xj , yj ) could be explained as the probability thatevent j is triggered by the ith event, ξ(tj , xj , yj ; ti , xi , yi ),when j > i,
(20)
ρij =
λ(tj , xj , yj )0,otherwise,
Moreover, the probability of the event j being a triggered event is
X
ρj =
ρij ,

(21)

i

and the probability that the jth event belongs to the background is
ϕj = 1 − ρ j =µ(xj , yj ).
λ(tj , xj , yj )(22)

If we select each event j with probabilities ρij , ρj or ϕj , we can form a newprocess being the process triggered by the ith event, the clustering process or the
background process, respectively. More details and an example are given in Section
4.
3.3 Estimating background rate
Given a dataset of earthquake occurrence times, locations and magnitudes in theobservation period T , the spatial distribution of the total seismicity rate is usuallyestimated by some nonparametric method, such as splines (Ogata 1998), kernel
functions (Zhuang et al. 2002; Helmstetter et al. 2007), grid averaging (Tsukakoshi
and Shimazaki 2006), or tessellation (Ogata 2004). In this paper, we consider thevariable (adaptive) kernel estimate following Zhuang et al. (2002). A similar adaptive
power law kernel was used by Werner et al. (2011). This adaptive approach issimple and tackles a serious disadvantage of the simple kernel estimate with a fixedbandwidth: for a spatially clustered point dataset, a small bandwidth gives a noisy
or variable estimate for the sparsely populated area, while a large bandwidth mixes
up the boundaries between the densely populated areas and the sparsely populatedareas. Therefore, instead of the kernel estimateNm̂1 (x, y) =1X
kd (x − xj , y − yj )
T j=1

(23)

Basic models of seismicity: Spatiotemporal models

9where kd (x, y) denotes the Gaussian kernel function
kd (x, y) =1 − x2 +y2 2
e 2d
2πd2

with a fixed bandwidth d, we adopt
N1Xm̂1 (x, y) =kd (x − xj , y − yj ),T j=1 j(24)

where dj represents the varying bandwidth calculated for each event j in the following way: Given a suitable integer np , find the smallest disk centered at the locationof the j-th event which includes at least np other earthquakes and is larger thansome minimum value (e.g., a distance within 0.02 degree, which is of the order ofthe location error, see also CORSSA Article in Theme IV) and let this minimum
radius be dj (e.g., Silverman 1986, Chapter 5).
Similar to estimating total seismicity, there have been many approaches for assessing background seismicity rate: (i) proportional to total seismicity rate of all
events or only of the big events in the catalog (Musmeci and Vere-Jones 1986; Console et al. 2003); (ii) using a declustering method to decluster the catalog and use thetotal rate in the declustered catalog as background rate (Ogata 1998); (iii) weighting each event by the probability that it is a background event (Zhuang et al. 2002,
2004); and (iv) Ogata (2004) introduced a Bayesian smoothness prior on tessellation
grids to estimate the spatial variation of the background and the model parametersat the same time. In this study, we use the third method because of its simplicity
together with unbiased estimates of sufficient precision.
Using the method described in the last paragraph of Section 3.2, once a background process is obtained, we can estimate the background intensity by applying asmoothing technique to the background catalog. Rather than repeating the thinning
procedure and the kernel estimation procedure many times to get an average estimate of the background intensity, we can directly estimate the average by weighting
all the events with their corresponding background probabilities, i.e.,µ̂(x, y) =1X
ϕi kdi (x − xi , y − yi ),
T i

(25)

where i runs over all of the events in the whole process, T is the length of process
duration, and kd is the Gaussian kernel function with a bandwidth d. The variablebandwidth di is defined in the same way as in (24).

10

www.corssa.org

3.4 Iterative algorithm
Sections 3.1 to 3.3 give us the method for estimating any one of the following three
things when the other two are known: the background rate µ, the model parameters
θ, and the background probabilities ϕi . In most cases, we have none of them but
a catalog of earthquakes. We can use the following algorithm to estimate themsimultaneously.
Algorithm A: stochastic declustering
A1. Given a fixed np and dmin , say 5 and 0.05 degree (roughly equivalent to 5.56 km
on Earth’s surface, which is close to the location error of earthquakes), calculate
the bandwidth dj for each event (tj , xj , Mj ), j = 1, 2, · · · , N .
A2. Set ` = 0, and u(0) (x, y) = 1.
A3. Using the maximum likelihood procedure (see, Ogata 1998), fit the model withconditional intensity functionX
λ(t, x, y) = νu(`) (x, y) +κ(Mi )g(t − ti )f (x − xi , y − yi ; Mi )(26)
i: ti <tto the earthquake data, where κ, g and f are defined in (10), and ν is therelaxation coefficient, introduced to speed up the convergence of the iterations.
A4. Calculate ρji , ρi and ϕi for each j < i and i = 1, 2, · · · , N .
A5. Calculate µ(x, y) from (25) and record it as u(`+1) (x, y).
A6. If max |u(`+1) (x, y) − u(`) (x, y)| > ε, where ε is a given small positive number(i.e., on the order of 10−6 ), then set ` = ` + 1 and go to step A3; otherwise, takeνu(`+1) (x, y) as the background rate and also output ρij , ρi and ϕi .
Example 1 Fitting ETAS model to the central Japan region
We use the earthquake data from the central-western Japan region from the
hypocenter catalog compiled by the Japan Meteorological Agency (JMA; see Fig.
1). We selected data for the period 1926-1995 in the rectangular area 34 ∼ 39◦ N
and 131 ∼ 140◦ E with magnitudes MJ ≥ 4.0 and depths ≤ 100km.
The space-time ETAS model equipped with the spatial response (18) is fit to theabove dataset. The MLEs for this region are Â = 0.020 event/(deg 2 ·day), α̂ = 1.135,ĉ = 0.040 days, p̂ = 1.15, and dˆ = .0010/deg 2 . Figures 1 (a)-(d) show estimates ofthe spatial variations of total intensity, clustering intensity, background intensity
and relative clustering effect, respectively. We can see that the clustering effects
arising from typical aftershocks are eliminated in the estimates of the background
intensity in Figure 1(d). In contrast, the background rates are almost the same as
the total intensities in the Wakayama and Ibaraki regions, which show almost purespatial occurrence of nonclustered events.

Basic models of seismicity: Spatiotemporal models
(a)

11

(b)
1000

40

100

39

1000

40

100

39

10
38

10
38

1
37

1
37

0.1
36

0.1
36

0.01
35

0.01
35

10−3

34

33
132

134

136

138

10−3

10−4

34

10−5

33

140

10−4

10−5
132

134

136

(c)

138

140

(d)
1000

40

100

39

40

1.0

39
0.8

10
38

38

1

0.6

37

37
0.1

36

36
0.01

0.4
35

35
10−3

34

33
132

134

136

138

140

10−4

34

10−5

33

0.2

132

134

136

138

140

≥ 4 earthquakes in central Japan (unit: events/(deg2 · year): (a) total, (b)Fig. 1 Estimated seismicity rate of Mj
clustered, and (c) background seismicity. (d) Clustering coefficient (ratio of clustered to total seismicity rate).

4 Stochastic reconstruction
Stochastic reconstruction was introduced to help classify earthquake events in a
catalog into different family trees. In testing hypotheses associated with earthquake
clustering features, such classification tackles the difficulties caused by the complicated mixture of the background seismicity and different earthquake clusters in both
space and time. This algorithm is outlined below.
Algorithm B: Stochastic classification of earthquake clusters [Zhuang et al. (2002),
Zhuang et al. (2004)].
B1. Calculate ϕi and ρji by using Algorithm A, where i = 1, 2, · · · , N and j =1, 2, · · · , i − 1. Here, N is the total number of events.
B2. For each event i, i = 1, 2, · · · , N , generate a random variable Ui uniformly
distributed on [0, 1].

12www.corssa.org

B3. For each i, let
Ii = min{k : ϕi +kXρji ≥ Ui and 0 ≤ k < i}.

j=1

If Ii = 0, then select i as a background or initial event; else, set the ith event tobe a direct offspring of the Ii th event.
We can repeat Algorithm B many times to get different stochastic versions ofseparations of the earthquake clusters. The non-uniqueness of such realizations illustrates the uncertainty in determining earthquake clusters, and thus repetition
can help us to evaluate the significance of some properties of seismicity clustering
patterns. However, we can also implement these tests by working with the probabilities ϕj and ρij directly. In the coming sections, we will show how to use theseprobabilities to reconstruct the characteristics associated with earthquake clusteringfeatures, using the Japanese JMA catalog as an example. The same reconstruction
procedures will be also applied to a simulated catalog for comparison.
Example 2 Stochastic declustering for the central Japan region Figure 2
shows a stochastically declustered version of the JMA catalog and illustrates thethinning procedure and the stochastic declustering algorithm, using the same data
as in Example 1.
Example 3 Location distributions of earthquake clusters
In the early studies of the space-time ETAS model, both scaled Gaussian spatial
response function (Model 2, Eq. 15) and scaled power-law spatial response (Model
4, Eq. 17) were adopted by Ogata (1998); Zhuang et al. (2002, 2004) as the spatial
response function. From the viewpoint of hazard mitigation, it is an importantquestion which of the above functions is appropriate for modeling seismicity, i.e.,
whether the aftershock activity decays quickly over short spatial ranges as in Model2, or decays in long distances as in Model 4. Ogata (1998) found that the second one
fit the data better than the first one by using the AIC model selection procedures.
Here we obtain the same conclusion from simple histogram techniques.
Define the transformed distance between a triggered event j and its direct ancestor, assumed i, ass(xj − xi )2 + (yj − yi )2
rij =
.
(27)
D2 exp[α(mi − m0 )]From (15) and (17), rij has a density function of2fR (r) = 2re−r , r ≥ 0;(28)

Basic models of seismicity: Spatiotemporal models

13

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 2 A realization of stochastic declustering for the central Japan region (MJ ≥ 4). (a), (b), and (c) are mapsof epicenter locations of all the earthquakes, the background events and the clustering events, respectively. (d), (e),and (f) are space-time plots of occurrence times against latitudes for all the earthquakes, the background eventsand the clustering events, respectively.

and2r(q − 1), r ≥ 0;(29)
(1 + r2 )q
for Model 2 and Model 4, respectively. The distribution with a density of (28)
is called a Rayleigh distribution. On the other hand, fR (r) can be reconstructed
throughPi,j ρij H(∆r/2 − |rij − r|)P,(30)fˆR (r) =∆r i,j ρijfR (r) =where ∆r is the bin width and H denotes the Heaviside function. The comparison
between fˆR and fR for the two models are shown in Figure 3. It can be seen that,when Model 2 is used, the reconstructed probability density of the transformed
distances between the ancestors and the direct offspring is quite different from thetheoretical one. When Model 4 is used, the reconstructed probability density isvery close to the theoretical one. These results confirm that the aftershocks show
a long range rather than a short range decay in space (Ogata 1998; Console and14www.corssa.org

0.8
0.6
0.0

0.2

0.4

p.d.f

0.6
0.4
0.0

0.2

p.d.f

0.8

1.0

MODEL 2

1.0

MODEL 1

0

2

4

6

8

10

Normalized distance

0

2

4

6

8

10

Normalized distanceFig. 3 Reconstruction results for the probability density functions (p.d.f.s) of the standardized triggering distancesfˆR (r) in equation (30) (circles) by using Model 2 (left panel) and Model 4 (right panel). The theoretical curves offR (r) in equations (28) and (29) are plotted in solid lines.

Murru 2001; Console et al. 2003). These results also imply the robustness of thereconstruction method because one gets a reconstructed probability density function
which is very close to the corresponding function in Model 4, even if an improper
model like Model 2 is employed.

5 Further reading on the ETAS model
1. Criticality. The ETAS model is a branching process, and thus its stability isalso related to the Rconcept of criticality. For the ETAS model, the criticality
parameter is ρ = κ(m) s(m) dm. When ρ < 1 the whole process is stable,
stationary and ergodic; if ρ ≥ 1, the population of the process in unit timeincreases monotonically with time and finally explodes. Detailed discussions canbe found in Helmstetter and Sornette (2002), Zhuang (2003), and Zhuang andOgata (2006).
2. Relation to foreshocks. It is an interesting question whether foreshocks can be explained by a clustering model for aftershocks, i.e., whether foreshocks are mainshocks whose aftershocks happen to be bigger. Studies by Helmstetter and Sornette (2003a), Helmstetter et al. (2003), Felzer et al. (2004), Zhuang and Ogata(2006), Zhuang et al. (2008), Christophersen and Smith (2008), Brodsky (2011)
and Marzocchi and Zhuang (2011) suggested that the probability of foreshock
phenomena does not exceed the probability of foreshock expected by an ETASlike generic clustering model.
3. Relation to Båth’s law. Båth’s law (Båth 1965) states that the median of themagnitude difference between a mainshock and its largest aftershock is aroundBasic models of seismicity: Spatiotemporal models

15

1.2. Utsu (1957) expressed the median of this difference as a linear function of themainshock magnitude. Helmstetter and Sornette (2003b); Saichev and Sornette(2005b); Zhuang and Ogata (2006); Vere-Jones and Zhuang (2008) showed in different ways that Båth’s law can be explained as a double exponential distribution
generated by the ETAS model.
4. HIST-ETAS model. Ogata et al. (2003) and Ogata (2004) developed a hierarchical space-time ETAS model to describe the difference of seismicity clustering
structures at different locations. He assumed that each parameter in the model
is a function of location, which has a smoothness prior. He developed powerful
Bayesian tools with penalized likelihoods to estimate the changes of clusteringcharacteristics in the form of spatial variation of model parameters.
5. Influence of small earthquakes. In the formulation of the ETAS model, the triggering effect from events smaller than the cut-off magnitude is ignored. The problem
was first studied by Sornette and Werner (2005a,2005b). Ignoring the triggering
effect may influence the conclusion of the study. Here we refer to Werner (2007);
Zhuang et al. (2008); Wang et al. (2010) for details.
6. Impact of anisotropy of earthquake clusters. Aftershock clusters of larger mainshocks are usually not isotropic around the epicenter of the mainshock but aligned
along the mainshock rupture. To account for anisotropy, Ogata (1998), Kagan andJackson (1995), Helmstetter et al. (2006) and Werner et al. (2011) used ellipticalaftershock distributions. Hainzl et al. (2008) demonstrated by means of ETAS
simulations that ignoring this fact can lead to a biased parameter estimation, inparticular, to an underestimation of the α-parameter.
7. Self-similar ETAS models. Vere-Jones (2005) developed a self-similar ETAS model
to avoid the problems caused by the cut-off magnitude in the ETAS model and enable fully self-similar features. Because of the missing data problem of immediate
aftershocks and small events, however, this model until now remains in theoretical developments (see also, Saichev and Sornette 2005a), and has not yet been
applied to real seismicity data. Another development to introduce self-similarity
is the branching aftershock sequence (BASS) model (e.g., Turcotte et al. 2007;
Holliday et al. 2008).
8. Earthquake forecasts. The ETAS model has been used for short-term (one-day)or real-time earthquake probability forecasts. There have been many versions,e.g., Vere-Jones (1998); Helmstetter et al. (2006); Werner et al. (2011); Zhuang(2010); Console et al. (2006, 2007, 2008, 2010); Marzocchi and Lombardi (2009);
Woessner et al. (2011).

16

www.corssa.org

6 The EEPAS model and the double branching model
The EEPAS (Every Earthquake is a Precursor According to Scale) was developedby Evison and Rhoades (2004) and Rhoades and Evison (2004) in the study offoreshock swarms. Its conditional intensity,Xλ(t, x, y, m) = µλ0 (t, x, y, m) +wi η(mi ) r(m|mi ) f (t − ti ) g(x − xi , y − yi | mi )i:ti <t(31)
has a similar form as the ETAS model but is very different in detail. The conditional
intensity function is not applied to the whole catalog but only for the events above
a second higher threshold, while the summation is taken over all the events above a
primary lower threshold. That is to say, this model could be regarded as a self- andmutually exciting process: the target process is of the events above the second higher
threshold, and the external process is of the events between the primary and thesecondary thresholds. In addition, f here takes the form of the lognormal density,and g and r take the form of normal densities, i.e.,f (u|mi ) =

uσT

1
√


(log u − aT − bT mi )2
exp −
,
2σT2
2π


(32)



1
x2 + y 2g(x, y|mi ) =exp − 2 mi b ,2πσA2 10mi bA
2σA 10 A(33)(m − aM − bM mi )2
exp −
,
2
2σM
2π

(34)

r(m|mi ) =

σM

1
√



where aT , bT , σT , bA , σA , aM , bM , and σM are model parameters to be estimated. The
weights wi are usually set to 1, but it has also been suggested to set wi to theprobability that the ith event is a background event, which could be obtained froman initial stochastic declustering using the ETAS model.
This model seems to produce high average probability gains or entropy scores formoderate term forecasts of moderate to large events. Here we refer to Rhoades andEvison (2004) for detailed applications.
The double branching model was proposed by Marzocchi and Lombardi (2008).
In this model, the triggering between earthquakes is divided into 2 steps: one is theshort-term clustering due to the elastic response of the upper layers of the Earth,which can be well represented by the ETAS model; and, the other is the interaction
among events over longer time-space scales. The conditional intensity of this modelBasic models of seismicity: Spatiotemporal models

17

takes the form,
λ(t, x, y) = µ2 (x, y) +X  K1 eα1 (Mi −M0 )i:ti <t+X
i:ti <t

K2 e

C1
2
p
(t − ti + c) (ri + d21 )q1

α2 (Mi −M0 ) −(t−ti )/τ

e



C2
2
(ri + d22 )q2(35)

where K1 , K2 , α1 , α2 , c, p, τ, d1 , RR
d2 , q1 and q2 are model parameters,
C2 are norRR C21 and
2
2 −q1
2 −q2
malizing constant such that C1 (ri +d1 ) dx dy = 1 and C2 (ri +d2 ) dx dy =1. Marzocchi and Lombardi (2008) used a two-step approach to estimate the model:
they firstly carried out Algorithm A to obtain the background probabilities by using the stochastic declustering method, and then fit a self-exciting model with thesecond summation term in (35) to the catalog of background events. They foundthat the background seismicity, obtained by removing short-term clustering , ischaracterized by a second-step long-term clustering with a characteristic time thatis compatible with post-seismic relaxation.
7 SummaryThis article presents an overview of the spatiotemporal models for seismicity, with
an emphasis on some powerful techniques associated with the ETAS model. As
can be seen in the article, the most important things are not what the model can
describe, but rather the phenomena that the model cannot describe. That is to say,one of the main outcomes of model building is determining aspects of the observed
process that the model cannot describe. We introduce the basic techniques of how
to decompose the earthquake clusters with uncertainty in the form of probabilities,
and how to use these probabilities to find the difference between real data and themodel. In summary, even though they are more difficult to implement and to apply
to earthquake data, spatiotemporal models are much more powerful in analyzing
seismicity than their simple temporal-only counterparts.
Acknowledgements The authors thank two anonymous reviewers for their constructive comments and Jeremy D. Zechar for editing this manuscript.
ReferencesBåth, M. (1965), Lateral inhomogeneities in the upper mantle, Tectonophysics, 2, 483–514. 14
Brodsky, E. E. (2011), The spatial density of foreshocks, Geophysical Research Letters, 38, L10,305, doi:10.1029/
2011GL047253. 14
Christophersen, A., and E. G. C. Smith (2008), Foreshock rates from aftershock abundance, Bulletin of the Seismological Society of Ameerica, 98 (5), 2133–2148, doi:10.1785/0120060143. 14

18

www.corssa.org

Console, R., and M. Murru (2001), A simple and testable model for earthquake clustering, Journal of GeophysicalResearch, 106, B5, 8699–8711. 13
Console, R., M. Murru, and A. M. Lombardi (2003), Refining earthquake clustering models, Journal of GeophysicalResearch, 108 (B10), 2468, doi:10.1029/2002JB002130. 6, 9, 14
Console, R., D. A. Rhoades, M. Murru, F. F. Evison, E. E. Papadimitriou, and V. Karakostas (2006), Comparative
performance of time-invariant, long-range and short-range forecasting models on the earthquake catalogue of
Greece, Journal of Geophysical Research, 111, B09,304, doi:10.1029/2005JB0044113. 15
Console, R., M. Murru, F. Catalli, and G. Falcone (2007), Real time forecasts through an earthquake clusteringmodel constrained by the rate-and-state constitutive law: comparison with a purely stochastic ETAS model,
Seismological Research Letters, 78 (1), 49–56. 3, 15
Console, R., M. Murru, and G. Falcone (2008), Real-time earthquake forecasting in Italy, Tectonophysics, 470 (3-4),
214–223. 15
Console, R., M. Murru, and G. Falcone (2010), Probability gains of an epidemic-type aftershock sequence modelin retrospective forecasting of m ≥ 5 earthquake in Italy, Journal of Seismology, 14 (1), 9–26, doi:10.1007/
s10950-009-9161-3. 15
Daley, D. D., and D. Vere-Jones (2003), An Introduction to Theory of Point Processes – Volume 1: Elementrary
Theory and Methods (2nd Edition), Springer, New York, NY. 7
Dieterich, J. (1994), A constitutive law for rate of earthquake production and its application to earthquake clustering,
Journal of Geophysical Research, 99 (B2), 2601–2618. 3
Evison, F. F., and D. A. Rhoades (2004), Demarcation and scaling of long-term seismogenesis, Pure and AppliedGeophysics, 161, 21–45, 10.1007/s00024-003-2435-8. 16
Felzer, K. R., R. E. Abercrombie, and G. Ekstrom (2004), A common origin for aftershocks, foreshocks, and multiplets, Bulletin of the Seismological Society of America, 94 (1), 88–98, doi:10.1785/0120030069. 14
Hainzl, S., A. Christophersen, and B. Enescu (2008), Impact of earthquake rupture extensions on parameter estimations of point-process models, Bulletin of the Seismological Society of Ameerica, 98 (4), 2066–2072, doi:
10.1785/0120070256. 15
Helmstetter, A., and D. Sornette (2002), Sub-critical and supercritical regimes in epidemic models of earthquake
aftershocks, Journal of Geophysical Research, 107 (B10), 2237. 14
Helmstetter, A., and D. Sornette (2003a), Foreshocks explained by cascades of triggered seismicity, Journal ofGeophysical Research, 108 (B10), 2457, doi:10.1029/2003JB002409. 14
Helmstetter, A., and D. Sornette (2003b), Båth’s law derived from the Gutenberg-Richter law and from aftershockproperties, Geophysical Research Letters, 30 (20), 2069, doi:10.1029/2003GL018186. 15
Helmstetter, A., D. Sornette, and J.-R. Grasso (2003), Mainshocks are aftershocks of conditional foreshocks: How
do foreshock statistical properties emerge from aftershock laws?, Journal of Geophysical Research, 108, 2046,
doi:10.1029/2002JB001991. 14
Helmstetter, A., Y. Y. Kagan, and D. D. Jackson (2006), Comparison of short-term and time-independent earthquake
forecast models for Southern California, Bulletin of the Seismological Society of Ameerica, 96 (1), 90–106, doi:
10.1785/0120050067. 15
Helmstetter, A., Y. Y. Kagan, and D. D. Jackson (2007), High-resolution time-independent forecast for m ≥ 5earthquakes in California, Seismogological Research Letters, 78, 59–67. 8
Holliday, J. R., D. L. Turcotte, and J. B. Rundle (2008), Self-similar branching of aftershock sequences, Physica A:
Statistical Mechanics and its Applications, 387 (4), 933 – 943, doi:DOI:10.1016/j.physa.2007.09.045. 15
Iwata, T. (2010), Effect of secondary aftershocks on aftershock decay in the rate- and state-friction model, Geophysical Research Abstracts, 12 (1), EGU2010–7920–2, EGU General Assembly. 3
Kagan, Y. (1991), Likelihood analysis of earthquake catalogues, Journal of Geophysical Research, 106, Ser. B7,
135–148. 4
Kagan, Y. Y., and D. D. Jackson (1995), New seismic gap hypothesis: Five years after, J. Geophys. Res., 100 (B3),3943–3959. 15
Lewis, P. A. W., and E. Shedler (1979), Simulation of non-nomogeneous Poisson processes by thinning, Naval
Research Logistics Quarterly, 26, 403–413. 7
Marzocchi, W., and A. M. Lombardi (2008), A double branching model for earthquake occurrence, Journal ofGeophysical Research, 113, B08,317, doi:10.1029/2007JB005472. 3, 16, 17Marzocchi, W., and A. M. Lombardi (2009), Real-time forecasting following a damaging earthquake, GeophysicalResearch Letters, 36, L21,302, doi:10.1029/2009GL040233. 15
Marzocchi, W., and J. Zhuang (2011), Statistics between mainshocks and foreshocks in Italy and Southern California,
Geophysical Research Letters, 38, L09,310, doi:10.1029/2011GL047165. 14

Basic models of seismicity: Spatiotemporal models

19

Musmeci, F., and D. Vere-Jones (1986), A variable-grid algorithm for smoothing clustered data, Biometrics, 42,
483–494. 9
Ogata, Y. (1981), On Lewis’ simulation method for point processes, IEEE Transactions on Information Theory,
IT-27 (1), 23–31. 7
Ogata, Y. (1988), Statistical models for earthquake occurrences and residual analysis for point processes, Journalof the American Statistical Association, 83, 9 – 27. 3
Ogata, Y. (1998), Space-time point-process models for earthquake occurrences, Annals of the Institute of Statistical
Mathematics, 50, 379–402. 4, 6, 7, 8, 9, 10, 12, 13, 15
Ogata, Y. (2004), Space-time model for regional seismicity and detection of crustal stress changes, Journal ofGeophysical Research, 109 (B3), B03,308, doi:10.1029/2003JB002621. 8, 9, 15
Ogata, Y., and J. Zhuang (2006), Space–time ETAS models and an improved extension, Tectonophysics, 413 (1-2),
13–23. 7
Ogata, Y., K. Katsura, and M. Tanemura (2003), Modelling heterogeneous space-time occurrences of earthquakes
and its residual analysis, Journal of the Royal Statistical Society: Series C (Applied Statistics), 52, 499–509(11),doi:doi:10.1111/1467-9876.00420. 15
Rathbun, S. L. (1993), Modeling marked spatio-temporal point patterns, Bulletin of the International Statistical
Institute, 55, Book 2, 379–396. 4, 6
Rhoades, D. A., and F. F. Evison (2004), Long-range earthquake forecasting with every earthquake a precursor
according to scale, Pure and Applied Geophysics, 161, 47–72, 10.1007/s00024-003-2434-9. 3, 16
Saichev, A., and D. Sornette (2005a), Vere-Jones’ self-similar branching model, Phys. Rev. E, 72 (5), 056,122, doi:
10.1103/PhysRevE.72.056122. 15
Saichev, A., and D. Sornette (2005b), Distribution of the largest aftershocks in branching models of triggered
seismicity: Theory of the universal Båth law, Physical Review, E, 71 (5), 056,127, doi:10.1103/PhysRevE.71.
056127. 15
Silverman, B. W. (1986), Density estimation for statistics and data analysis, Chapman and Hall,, London, New
York. 9Sornette, D., and M. J. Werner (2005a), Apparent clustering and apparent background earthquakes biased byundetected seismicity, J. Geophys. Res. Solid Earth, 110 (B09303), doi:10.1029/2005JB003621. 15
Sornette, D., and M. J. Werner (2005b), Constraints on the size of the smallest triggering earthquake from theepidemic-type aftershock sequence model, Båth’s law, and observed aftershock sequences, Journal of GeophysicalResearch, 110, B08,304, doi:10.1029/2004JB003535. 15
Tsukakoshi, Y., and K. Shimazaki (2006), Temporal behavior of the background seismicity rate in central Japan,
1998 to mid-2003, Tectonophysics, 417, 155–168. 8
Turcotte, D., J. Holliday, and J. Rundle (2007), BASS, an alternative to ETAS, Geophysical Research Letters,
34 (12), cited By (since 1996) 8. 15
Utsu, T. (1957), Magnitude of earthquakes and occurrence of their aftershocks, Zisin (J. Seismol. Soc. Jap.), 10,
35–45 (in Japanese). 15
Veen, A., and F. P. Schoenberg (2008), Estimation of space-time branching process models in seismology using an EM-type algorithm, Journal of the American Statistical Association, 103 (482), 614–624, doi:10.1198/
016214508000000148. 7
Vere-Jones, D. (1998), Probability and information gain for earthquake forecasting, Computational Seismology, 30,
248–263. 15
Vere-Jones, D. (2005), A class of self-similar random measure, Advances in Applied Probability, 37 (4), pp. 908–914.
15
Vere-Jones, D., and J. Zhuang (2008), Distribution of the largest event in the critical epidemic-type aftershocksequence model, Phys. Rev. E, 78 (4), 047,102, doi:10.1103/PhysRevE.78.047102. 15
Wang, Q., D. D. Jackson, and J. Zhuang (2010), Missing links in earthquake clustering models, Geophysical Research
Letters, 37, L21,307. 7, 15
Werner, M. J. (2007), On the fluctuations of seismicity and uncertainties in earthquake catalogs: Implications andmethods for hypothesis testing, Ph.D. thesis, Univ. of Calif., Los Angeles. 15
Werner, M. J., A. Helmstetter, D. D. Jackson, and Y. Y. Kagan (2011), High-Resolution Long-Term and Short-Term
Earthquake Forecasts for California, Bull. Seism. Soc. Am., 101 (4), 1630–1648, doi:10.1785/0120090340. 8, 15
Woessner, J., S. Hainzl, W. Marzocchi, M. J. Werner, A. M. Lombardi, F. Catalli, B. Enescu, M. Cocco, M. C.
Gerstenberger, and S. Wiemer (2011), A retrospective comparative forecast test on the 1992 Landers sequence,
Journal of Geophysical Research, 116, B05,305, doi:10.1029/2010JB007846. 15

20

www.corssa.org

Zhuang, J. (2003), Some applications of point processes in seismicity modelling and prediction, Ph.D. thesis, theGraduate Unversity for Advanced Studies, Kanagawa, Japan. 14
Zhuang, J. (2006), Second-order residual analysis of spatiotemporal point processes and applications in model
evaluation, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68 (4), 635–653, doi:
10.1111/j.1467-9868.2006.00559.x. 7
Zhuang, J. (2010), Next-day earthquake forecasts by using the etas model, submitted to Earth, Planet, and Space,in press. 15
Zhuang, J., and Y. Ogata (2006), Properties of the probability distribution associated with the largest event in an
earthquake cluster and their implications to foreshocks, Physical Review, E, 73, 046,134, doi:10.1103/PhysRevE.
73.046134. 14, 15
Zhuang, J., Y. Ogata, and D. Vere-Jones (2002), Stochastic declustering of space-time earthquake occurrences,
Journal of the American Statistical Association, 97 (3), 369–380. 6, 8, 9, 11, 12
Zhuang, J., Y. Ogata, and D. Vere-Jones (2004), Analyzing earthquake clustering features by using stochastic
reconstruction, Journal of Geophysical Research, B5 (3), B05,301, doi:10.1029/2003JB002879. 6, 7, 9, 11, 12
Zhuang, J., C.-P. Chang, Y. Ogata, and Y.-I. Chen (2005), A study on the background and clustering seismicity
in the Taiwan region by using a point process model, Journal of Geophysical Research, 110, B05S13, doi:
10.1029/2004JB003157. 7
Zhuang, J., A. Christophosen, M. K. Savage, D. Vere-Jones, Y. Ogata, and D. D. Jackson (2008), Differences between spontaneous and triggered earthquakes: their influences on foreshock probabilities, Journal of GeophysicalResearch, 113, B11,302, doi:10.1029/2008JB005579. 14, 15