tŚĂƚŝƐƚŚĞŚĂǌĂƌĚĨŽƌƚŚŝƐƐŝƚĞŝŶZĞŶŽ͍͞^ĞŝƐŵŝĐŚĂǌĂƌĚĂŶĂůǇƐĞƐŝŶƚŚƌĞĞĚĞĐĂĚĞƐŽĨƌĞƐĞĂƌĐŚǁŝƚŚ ƐƉĞĐŝĂůĞŵƉŚĂƐŝƐŽŶƚŚĞŝŵƉŽƌƚĂŶĐĞŽĨƋƵĂŶƚŝƚĂƚŝǀĞƐƚƌŽŶŐŵŽƚŝŽŶƉƌĞĚŝĐƚŝŽŶƐ͟ EƵŵĞƌŽƵƐŚŝƐƚŽƌŝĐĂůĞĂƌƚŚƋƵĂŬĞƐ͘ ĐƚŝǀĞĨĂƵůƚƐŶĞĂƌďǇ͘:ŽŚŶ'͘ŶĚĞƌƐŽŶ
WƌŽĨĞƐƐŽƌŽĨ'ĞŽƉŚǇƐŝĐƐ
hŶŝǀĞƌƐŝƚǇŽĨEĞǀĂĚĂ͕ZĞŶŽ'ƌĞĞŶ͗DхϰzĞůůŽǁ͗DхϱZĞĚ͗DхϲŽŶĐĞƉƚƵĂůƉƌŽĐĞƐƐ YƵĂŶƚŝĨǇƚŚĞŚĂǌĂƌĚĨƌŽŵĞĂĐŚƐŽƵƌĐĞƐĞƉĂƌĂƚĞůǇ͘ ĚĚƚŚĞŚĂǌĂƌĚĨƌŽŵĂůůƐĞƉĂƌĂƚĞƐŽƵƌĐĞƐƚŽŽďƚĂŝŶĂĨŝŶĂůƌĞƐƵůƚ͘

ϭ

Ϯ

−37−

,ĂǌĂƌĚƵƌǀĞEĂƚŝŽŶĂů^ĞŝƐŵŝĐ,ĂǌĂƌĚDĂƉ ƐƐƵŵƉƚŝŽŶƐ͗ ƚƌƵĞ͕ĚĂƚĂͲŐĞŶĞƌĂƚŝŶŐŚĂǌĂƌĚĐƵƌǀĞĞǆŝƐƚƐĨŽƌĞǀĞƌǇƐŝƚĞ͘ DŽĚĞůƐĂƌĞƌĞƋƵŝƌĞĚƚŽĞƐƚŝŵĂƚĞƚŚĞŚĂǌĂƌĚĐƵƌǀĞ͘ZĞŶŽ ŽŵŵƵŶŝƚǇĞĨĨŽƌƚ h^'^ĞŶĐŽƵƌĂŐĞƐĂŶĚƐĞĞŬƐŝŶƉƵƚĨƌŽŵƐĞŝƐŵŽůŽŐŝƐƚƐĂŶĚƵƐĞƌƐ
DĂƉŐĞŶĞƌĂƚĞĚďǇĐŽŶƚŽƵƌŝŶŐƚŚĞŐƌŽƵŶĚŵŽƚŝŽŶƐĨƌŽŵŚĂǌĂƌĚĐƵƌǀĞƐŐĞŶĞƌĂƚĞĚĨŽƌĂĚĞŶƐĞŐƌŝĚŽĨƉŽŝŶƚƐ͘'DWD6RXUFH86*6ZHEVLWH

ϯ^ŽƵƌĐĞ͗h^'^K&ZϮϬϭϰͲϭϬϵϭϰ^ǇŵďŽůŝĐĞƋƵĂƚŝŽŶƚŽďƵŝůĚĂŚĂǌĂƌĚĐƵƌǀĞ

KƵƚůŝŶĞ /ĚĞŶƚŝĨǇϭϱŬĞǇƉƌŽďůĞŵƐƚŚĂƚŶĞĞĚƚŽďĞƐŽůǀĞĚĨŽƌƐĞŝƐŵŝĐŚĂǌĂƌĚƐƚƵĚŝĞƐ͘ &ŽƌƐŽŵĞ͗ũƵƐƚĚĞƐĐƌŝďĞƚŚĞƉƌŽďůĞŵ͘ KƚŚĞƌƐ͗ĚŝƐĐƵƐƐŝŶŵŽƌĞĚĞƚĂŝů͘zсŐƌŽƵŶĚŵŽƚŝŽŶ
γ;zͿ = ƚŚĞŚĂǌĂƌĚĐƵƌǀĞ͕ƚŚĞƌĂƚĞŽĨĞǆĐĞĞĚŝŶŐz͘ŶсŶƵŵďĞƌŽĨĞĂƌƚŚƋƵĂŬĞƐ͕ŵĂŐŶŝƚƵĚĞD͕ĂƚůŽĐĂƚŝŽŶǆΦ сƉƌŽďĂďŝůŝƚǇŽĨƚŚĞĞĂƌƚŚƋƵĂŬĞĞǆĐĞĞĚŝŶŐz͘
ϱ

ϲ

−38−

^ŚŽƵůĚƚŚĞƐĞŝƐŵŝĐŝƚǇĐĂƚĂůŽŐďĞĚĞĐůƵƐƚĞƌĞĚ͍/ĨƐŽ͕ŚŽǁ͍^ĞŝƐŵŝĐŝƚǇDŽĚĞů
ĞĐůƵƐƚĞƌĞĚ ƐĞŝƐŵŝĐŝƚǇŵŽĚĞů WƵƌƉŽƐĞ͗ƚŽŵĂŬĞƚŚĞĐĂƚĂůŽŐWŽŝƐƐŽŶŝĂŶ͘ ^ŽŵĞŝŶĨŽƌŵĂƚŝŽŶŝƐŶŽƚƵƐĞĚ͘tŚĂƚĂƌĞƚŚĞƚƌĂĚĞŽĨĨƐ͍&ĂƵůƚƐ

ϳϴtŚĂƚŝƐƚŚĞŵĂŐŶŝƚƵĚĞʹ ĨƌĞƋƵĞŶĐǇĚŝƐƚƌŝďƵƚŝŽŶĨŽƌƚŚĞŶĞƚǁŽƌŬŽĨĨĂƵůƚƐ͍>ĨсϯϲŬŵŝƉсϱϬŽtĨсϭϵ͘ϲŬŵ^ĨсϮŵŵͬǇƌ ƐĐĂůŝŶŐŵŽĚĞůŝƐƵƐĞĚƚŽĞƐƚŝŵĂƚĞDǁ͘^ŚŽƵůĚŝƚ͗ hƐĞĨƵůůůĞŶŐƚŚ͍ hƐĞĂĨƌĂĐƚŝŽŶŽĨůĞŶŐƚŚ͍
 >ŝŵŝƚůĞŶŐƚŚďǇďĞŶĚƐ͕ƐƚĞƉƐ͕͙͍ hZ&ͲϯĂƉƉƌŽĂĐŚ͍ KƚŚĞƌhŶĐĞƌƚĂŝŶƚŝĞƐ͗ĚŝƉ͕ĚŽǁŶĚŝƉŐĞŽŵĞƚƌǇ͕ŵĂǆŝŵƵŵĚĞƉƚŚŽĨĨĂƵůƚŝŶŐ͘ǆĂŵƉůĞ͗DŽƵŶƚZŽƐĞ&ĂƵůƚŽŶĞ͗&Ƶůů&ĂƵůƚZƵƉƚƵƌĞĂƌƚŚƋƵĂŬĞƉƉƌŽĂĐŚ
^ĐĂůŝŶŐZĞůĂƚŝŽŶƐŚŝƉ͗>Ĩ Æ DǁtĞůůƐΘŽƉƉĞƌƐŵŝƚŚ;ϭϵϵϰͿ
Dǁсϲ͘ϴϮ
DŽŵĞŶƚсϮ͘ϭϴΎϭϬϮϲ ĚǇŶĞͲĐŵ^ůŝƉсDŽŵĞŶƚͬ;ƐŚĞĂƌŵŽĚƵůƵƐΎĨĂƵůƚĂƌĞĂͿсϭ͘ϭϲŵƋ͘ƌĂƚĞс;ƐůŝƉƌĂƚĞͬƐůŝƉƉĞƌĞǀĞŶƚͿсϭͬϱϴϬсϬ͘ϬϬϭϳĞǀĞŶƚƐͬǇĞĂƌ^ŽŶ;DͿĨŽƌƚŚŝƐĨĂƵůƚсŶ;ϲ͘ϴϮͿсϬ͘ϬϬϭϳǇƌ Ͳϭ

ϵ

ϭϬ

−39−

&Ƶůů&ĂƵůƚZƵƉƚƵƌĞĂƌƚŚƋƵĂŬĞƉƉƌŽĂĐŚ
>ĨсϯϲŬŵŝƉсϱϬŽtĨсϭϵ͘ϲŬŵ^ĨсϮŵŵͬǇƌ&ĂŝůƵƌĞƐŽĨ&Ƶůů&ĂƵůƚZƵƉƚƵƌĞDŽĚĞů
ϭ͘͞^ŵĂůů͟ĞĂƌƚŚƋƵĂŬĞƐͲ ƐŽŵĞŵĂǇŚĂǀĞďĞĞŶŽŶƚŚŝƐĨĂƵůƚ͘ϱ&ĂƵůƚƐÆ ůŝƐƚŽĨϱŵĂŐŶŝƚƵĚĞͲƌĂƚĞƉĂŝƌƐ

DŽƵŶƚZŽƐĞDǁсϲ͘ϴϮϬ͘ϬϬϭϳ

<ŝŶŐƐĂŶǇŽŶ

ϲ͘ϱϮϬ͘ϬϬϬϰϬ

ĂƌƐŽŶŝƚǇϲ͘ϰϴϬ͘ϬϬϬϮϭ

/ŶĚŝĂŶ,ŝůůϲ͘ϭϯϬ͘ϬϬϬϯϱ

ĂƌƐŽŶZĂŶŐĞ

ϳ͘ϬϴϬ͘ϬϬϭϳ

Ϯ͘ĂƐŬĞǇĞƚĂů͘;^^͕ϭϵϵϲͿĨŝŶĚŝŶŐƐĨŽƌϭϵϱϰ&ĂŝƌǀŝĞǁWĞĂŬĂŶĚŝǆŝĞsĂůůĞǇĞĂƌƚŚƋƵĂŬĞƐ͗^ĞǀĞƌĂůĨĂƵůƚƐǁĞƌĞŝŶǀŽůǀĞĚ͘ĂƐŬĞǇ͕^͘:͕͘ĞƚĂů͘;ϭϵϵϲͿ
Ƶůů͘^ĞŝƐŵ͘^ŽĐ͘ŵ͕͘ǀ͘ϴϲ͕Ɖ͘ϳϲϭͲ
ϳϴϳ͘

ϭϭ

/ŵƉůŝĐĂƚŝŽŶʹ ǁĞŶĞĞĚƚŽĐŽŶƐŝĚĞƌŵƵůƚŝƉůĞĨĂƵůƚƌƵƉƚƵƌĞƐŝŶƚŚĞƐĞŝƐŵŝĐŚĂǌĂƌĚŵŽĚĞůƐ͘ϭϮDŽĚŝĨŝĞĚhZ&ϯƉƉƌŽĂĐŚ

,ŽǁĚĞĞƉůǇĚŽůĂƌŐĞĐŽŶƚŝŶĞŶƚĂůĞĂƌƚŚƋƵĂŬĞƐƌƵƉƚƵƌĞ͍/ƐƚŚĞĚĞƉƚŚƉƌĞĚŝĐƚĂďůĞ͍ ^ƵďĚŝǀŝĚĞĞĂĐŚĨĂƵůƚŝŶƚŽƐƵďƐĞĐƚŝŽŶƐ͘ ƌĞĂƚĞĂƐĞƚŽĨ͞ƉŽƐƐŝďůĞ͟ƌƵƉƚƵƌĞƐƵƐŝŶŐĐŽŵďŝŶĂƚŝŽŶƐŽĨƐƵďƐĞĐƚŝŽŶƐ͘ &ŝŶĚƌĂƚĞƐĨŽƌĞĂĐŚƌƵƉƚƵƌĞʹ ďƵƚƚŚĞƐŽůƵƚŝŽŶŝƐŚŝŐŚůǇƵŶĚĞƌĚĞƚĞƌŵŝŶĞĚ͘ /ŵƉŽƌƚĂŶƚĨŽƌ͗ ^ŽƵƌĐĞƐĐĂůŝŶŐ;ĞƐƚŝŵĂƚĞDǁĨƌŽŵŽďƐĞƌǀĂďůĞĨĂƵůƚůĞŶŐƚŚͿ ĂůĂŶĐŝŶŐƐĞŝƐŵŝĐŵŽŵĞŶƚ ŽŵŝĐƌŽĞĂƌƚŚƋƵĂŬĞƐƌĞǀĞĂůƚŚĞĚĞƉƚŚŽĨƌƵƉƚƵƌĞ͍

ϭϯ

ϭϰ

−40−

/ƐŵŽŵĞŶƚďĂůĂŶĐŝŶŐďĂƐĞĚŽŶƐůŝƉͬƐƚƌĂŝŶƌĂƚĞƚŚĞďĞƐƚĂƉƉƌŽĂĐŚƚŽĞƐƚŝŵĂƚĞĞĂƌƚŚƋƵĂŬĞƌĂƚĞƐ͍ƌĞƚŚĞƌĞŽƚŚĞƌŵĞƚŚŽĚƐ͍ ŽĞƐŶ͛ƚǁŽƌŬĨŽƌůŽǁͲŚĂǌĂƌĚĂƌĞĂƐ͕ƐƵĐŚĂƐĞĂƐƚĞƌŶEŽƌƚŚŵĞƌŝĐĂ͘ tŚĂƚĐŽŶƚƌŽůƐƚŚĞƌĂƚĞŽĨĞĂƌƚŚƋƵĂŬĞƐŝŶ͞ƐƚĂďůĞĐŽŶƚŝŶĞŶƚĂůƌĞŐŝŽŶƐ͍͟ YƵĂŶƚŝĨǇŽĨĨͲĨĂƵůƚƐƚƌĂŝŶ͍ ƌĞƚŚĞƌĞĂůƚĞƌŶĂƚŝǀĞƐƚŽďĂĐŬŐƌŽƵŶĚƐĞŝƐŵŝĐŝƚǇƚŽƋƵĂŶƚŝĨǇƚŚĞƌŽůĞŽĨŵŝŶŽƌ͕ƵŶŵĂƉƉĞĚĨĂƵůƚƐ͍'ƌŽƵŶĚDŽƚŝŽŶDŽĚĞůƐ;'DWƐͿ 'ƌŽƵŶĚDŽƚŝŽŶWƌĞĚŝĐƚŝŽŶƋƵĂƚŝŽŶƐ

&DPSEHOO %R]RUJQLD 

3*$J ^ǇŶƚŚĞƚŝĐ^ĞŝƐŵŽŐƌĂŵƐ 
0 : 
0 : 
0 :  ZĞĨĞƌĞŶĐĞƐ͗ŶĚĞƌƐŽŶ;ϭϵϳϵͿ͕ŝƌĚΘ<ƌĞĞŵĞƌ ;ϮϬϭϱͿ

0 : 

 
 

 

 

U UXS

ϭϱ

ϭϲ

tŚĂƚŝƐƚŚĞŐƌŽƵŶĚŵŽƚŝŽŶǀĞƌǇŶĞĂƌDǁхϳĞĂƌƚŚƋƵĂŬĞƐŝŶƚŚĞĐŽŶƚŝŶĞŶƚĂůĐƌƵƐƚ͍E'tĞƐƚϮĂƚĂďĂƐĞ ĂƚĂ͍ &ĂƌŵŽƌĞƚŚĂŶϯϬǇĞĂƌƐĂŐŽ͕ďƵƚŶŽƚĞŶŽƵŐŚĨŽƌůĂƌŐĞD͕ƐŚŽƌƚZ͘
 &ŽƌĚŝƉͲƐůŝƉĞǀĞŶƚƐ͗ǁŚĂƚŚĂƉƉĞŶƐŽŶƚŚĞŚĂŶŐŝŶŐǁĂůůĂŶĚĨŽŽƚǁĂůů͍ ZŝƐŬŽĨĐŽŶƐŝĚĞƌŝŶŐĂŶǇĞĂƌƚŚƋƵĂŬĞǁĞĚŽƐĞĞĂƐ͞ƚǇƉŝĐĂů͟ dŚĞƐĞĞǀĞŶƚƐĂƌĞƌĂƌĞ͘ ^ƚĂƚƵƐŽĨŝŶƐƚƌƵŵĞŶƚĂƚŝŽŶ ǀĞƌǇŽŶĞƐŚŽƵůĚĚŽǁŚĂƚƚŚĞǇĐĂŶʹ ŵŝƐƐĞĚŽƉƉŽƌƚƵŶŝƚŝĞƐĂƌĞĐŽƐƚůǇ͘ WƵƌƐƵĞƉŽƐƐŝďůĞĨŽƌĞƐŚŽĐŬƐ͘

ϭϳ

ϭϴ

−41−

ŚĂŶĐĞƐŽĨŐĞƚƚŝŶŐĚĂƚĂZŝƐŬŽĨŽǀĞƌƌĞůŝĂŶĐĞŽŶĨĞǁĞĂƌƚŚƋƵĂŬĞƐ
 'DWƐŵĂǇĂƐƐƵŵĞŵŽƚŝŽŶƐƐĂƚƵƌĂƚĞŶĞĂƌĨĂƵůƚʹ ŵŽƚŝŽŶĚŽŵŝŶĂƚĞĚďǇŶĞĂƌďǇƐůŝƉ͘ ŽƵŶƚĞƌĞǆĂŵƉůĞŝŶϭϵϴϱDŝĐŚŽĂĐĂŶ͕DĞǆŝĐŽĞĂƌƚŚƋƵĂŬĞ;Dϴ͘ϭͿĂŶĚϮϬϭϭdŽŚŽŬƵĞĂƌƚŚƋƵĂŬĞ;Dϵ͘ϭͿ͘ ^D'ϯǁĂƐΕϭϬǆŵŽƌĞƉŽǁĞƌĨƵů͘^ƚƌŝŬĞͲ^ůŝƉ

Dхϳ

Dхϳ͘ϱ

Ϭ͘ϰ;ϭͬϮ͘ϱǇͿϬ͘ϭϱ;ϭͬϳǇͿZĞǀĞƌƐĞ

Ϭ͘Ϭϳ;ϭͬϭϱǇͿϬ͘Ϭϰϱ;ϭͬϮϮǇͿ

EŽƌŵĂů

Ϭ͘Ϭϰϱ;ϭͬϮϮǇͿEŽŶĞǇĞƚ

hƉĚĂƚĞĚĨƌŽŵŶĚĞƌƐŽŶĞƚĂů͘;ϮϬϭϳ͕^^ͿϭϵϮϬs^ϯϬ ŝƐĂƉƌŽǆǇ
/ŶE'tĞƐƚϮ&ůĂƚĨŝůĞ

,ŽǁƐŚŽƵůĚƐŝƚĞƌĞƐƉŽŶƐĞďĞƌĞƉƌĞƐĞŶƚĞĚŝŶ
'DWƐ͍ ϯϬǇĞĂƌƐĂŐŽǁĞǁĞƌĞƵƐŝŶŐƋƵĂůŝƚĂƚŝǀĞƚĞƌŵƐƚŽĚĞƐĐƌŝďĞƚŚĞƐŝƚĞƉƌŽƉĞƌƚŝĞƐ͘ ǆƚĞŶƐŝǀĞŵĞĂƐƵƌĞŵĞŶƚƐƐŝŶĐĞƚŚĞŶ͘ ƵƌƌĞŶƚůǇ͕s^ϯϬ ŝƐĂĐŽŵŵŽŶƉƌŽǆǇĨŽƌƐŝƚĞƌĞƐƉŽŶƐĞŝŶ
'DWƐ͘ ĂŶǁĞĚŽďĞƚƚĞƌ͍,ĂƌĚƐŝƚĞƐŽŶƚŽƉ^ŽĨƚƐŝƚĞƐŽŶƚŽƉ

ϮϮ

Ϯϭ

−42−
$SULO)XNXVKLPD+DPDGRULHDUWKTXDNH0:

$SULO)XNXVKLPD+DPDGRUL HDUWKTXDNH0:
,PDJHVIURP0DUX\DPDHWDO *6-

ĂďŝŶĂĚũĂĐĞŶƚƚŽĨĂƵůƚƚŚĂƚƌƵƉƚƵƌĞĚŝŶƚŚĞϭϵϱϰŝǆŝĞsĂůůĞǇĞĂƌƚŚƋƵĂŬĞ;Dtϲ͘ϵͿ͘dŚĞƉŚŽƚŽŐƌĂƉŚǁĂƐƚĂŬĞŶŽŶϮϬĞĐĞŵďĞƌϭϵϱϰďǇ<Ăƌů^ƚĞŝŶďƌƵŐŐĞ ;E/^͕hͿ͘0DUFK0:,PPHGLDWHDIWHUVKRFNV&ĂƵůƚĞǆƉŽƐƵƌĞĚƵƌŝŶŐƚŚĞϮϬϭϭ&ƵŬƵƐŚŝŵĂ,ĂŵĂĚŽƌŝ ĞĂƌƚŚƋƵĂŬĞ͘dŚĞƐĞƉŚŽƚŽŐƌĂƉŚƐǁĞƌĞƚĂŬĞŶďǇdĂĚĂƐŚŝDĂƌƵǇĂŵĂ͕'ĞŽůŽŐŝĐĂů^ƵƌǀĞǇŽĨ:ĂƉĂŶ͘

Ϯϯ

Ϯϰ

ď

Ă

ϯϲϬͲϳϲϬ

хϳϲϬ<ͲEdĂŶĚ<ŝ<ͲŶĞƚ/ŶĨůƵĞŶĐĞŽĨs^ϯϬ ŝƐŵŽƌĞĐŽŶƐƉŝĐƵŽƵƐ͘ĐĚϭϴϬͲϯϲϬ

фϭϴϬ&ƵŬƵƐŚŝŵĂͲ,ĂŵĂĚŽƌŝ
ϭϭƉƌŝů͕ϮϬϭϭ
Dǁϲ͘ϳ

Ϯϱ

ĂŶƐŝŐŵĂ;σͿďĞƌĞĚƵĐĞĚŝŶ'DWƐ͍ ƌŝƚŝĐĂůƌŽůĞŽĨƐŝŐŵĂĨŽƌW^,ƌĞƐƵůƚƐĂƚůŽǁƉƌŽďĂďŝůŝƚŝĞƐ͘ ƌŐŽĚŝĐĂƐƐƵŵƉƚŝŽŶ;ŶĚĞƌƐŽŶΘƌƵŶĞ͕ϭϵϵϵ͕^Z>Ϳ
 ^ŝŶŐůĞͲƐƚĂƚŝŽŶƐŝŐŵĂ

Ϯϲ

−43−
 σ ŝƐƚŚĞŵĞĂƐƵƌĞŽĨƚŚĞŵŝƐĨŝƚŽĨĚĂƚĂĨƌŽŵĂ'DW͘ &ŝŐƵƌĞŚĂƐĂƉƉůŝĞĚĂĚũƵƐƚŵĞŶƚƚŽƉĞĂŬǀĞůŽĐŝƚǇ;W'sͿĨƌŽŵ<ĂǁĂƐĞ
ĂŶĚDĂƚƐƵŽ;ϮϬϬϰĂ͕ďͿ͘ ^ĐĂƚƚĞƌƌĞŵĂŝŶƐ͘ ƌŐŽĚŝĐĂƐƐƵŵƉƚŝŽŶ͗ƵƐĞƐƚŚŝƐƐƉĂƚŝĂůƐĐĂƚƚĞƌĂƐĂŵĞĂƐƵƌĞŽĨƚŚĞƵŶĐĞƌƚĂŝŶƚǇĂƚĂƐŝŶŐůĞůŽĐĂƚŝŽŶŝŶƚŝŵĞ͘ /ŵƉůŝĐĂƚŝŽŶ͗ĨŽƌƐĞŝƐŵŝĐŚĂǌĂƌĚǁĞƐŚŽƵůĚĚĞƚĞƌŵŝŶĞσ ĨŽƌƚŚĞůŽĐĂƚŝŽŶŽĨƚŚĞƐƚƌƵĐƚƵƌĞ͘ϮϳWƌĞĐĂƌŝŽƵƐƌŽĐŬƐĂƚ>ŽǀĞũŽǇƵƚƚĞ͕͕ĂƌĞŝŶĐŽŶƐŝƐƚĞŶƚǁŝƚŚƚŚĞŚĂǌĂƌĚŵŽĚĞů͘/Ĩσ ŝƐƐŵĂůůĞƌ͕ƚŚĞǇĂƌĞĐŽŶƐŝƐƚĞŶƚ͘EŽƚĞ͗ĂůƐŽƵƐŝŶŐƉƌĞĐĂƌŝŽƵƐƌŽĐŬƐƚŽŚĞůƉĐĂůŝďƌĂƚĞƐǇŶƚŚĞƚŝĐƐ͘Ϯϴ,ŽǁĚŽǁĞĂĚũƵƐƚŐƌŽƵŶĚŵŽƚŝŽŶŵŽĚĞůƐĨŽƌĚŝĨĨĞƌĞŶƚƌĞŐŝŽŶĂůŐĞŽůŽŐǇΘƚĞĐƚŽŶŝĐƐ͍ǆĂŵƉůĞ͗&ƵŬƵƐŚŝŵĂ,ĂŵĂĚŽƌŝ ƌĞĐŽƌĚƐĞǆĐĞĞĚh^'DWƐ͘ŽŶƐŝƐƚĞŶƚĨŽƌW'ĂŶĚW'sĨŽƌŵŽƐƚĞǀĞŶƚƐ;DхϱͿŝŶƚŚĞƐĞƋƵĞŶĐĞ͘ ZĞŐŝŽŶĂůĚŝĨĨĞƌĞŶĐĞƐŝŶƚĞĐƚŽŶŝĐƐ͕ĐƌƵƐƚĂůƚŚŝĐŬŶĞƐƐ͕ƚĞŵƉĞƌĂƚƵƌĞ͕ƌŽĐŬƚǇƉĞ͘ ZĞƐƵůƚ͗'DWƐŽƌƐǇŶƚŚĞƚŝĐƐƐĞŝƐŵŽŐƌĂŵƐŵŝƐĐĂůŝďƌĂƚĞĚ͘

Ϯϵ

ϯϬ

−44−

tŚĂƚŝƐƚŚĞƐŚĂƉĞŽĨƚŚĞƐĞŝƐŵŝĐ͞ƐŽƵƌĐĞƐƉĞĐƚƌƵŵ͕͟;ĨͿ͍^ǇŶƚŚĞƚŝĐ^ĞŝƐŵŽŐƌĂŵƐ

 /ŶǀĞƌƐŝŽŶƐĨŽƌƐŽƵƌĐĞŵŽĚĞůƐĐŚĂƌĂĐƚĞƌŝǌĞ;ĨͿĂƚůŽǁĨƌĞƋƵĞŶĐŝĞƐ͘ /ƐŽŵĞŐĂͲƐƋƵĂƌĞĂŐŽŽĚĂǀĞƌĂŐĞĂƚŚŝŐŚĨƌĞƋƵĞŶĐŝĞƐ͍ ŽĞƐ;ĨͿĚĞƉĞŶĚŽŶĚĞƉƚŚ͍͙ƌŚĞŽůŽŐǇ͍ϯϭϯϮĐĐĞůĞƌŽŐƌĂŵĨƌŽŵDǁϲ͘ϳ&ƵŬƵƐŚŝŵĂͲ,ĂŵĂĚŽƌŝ ϭϭƉƌ͘ϮϬϭϭ ZŝŐŚƚ͗&ŽƵƌŝĞƌĂŵƉůŝƚƵĚĞŽĨĞĂƐƚĐŽŵƉŽŶĞŶƚ͕ƐƵƌĨĂĐĞΘĚŽǁŶŚŽůĞ
 ,ŝŐŚĨƌĞƋƵĞŶĐǇ͞ĨůĂƚ͟ƐƉĞĐƚƌƵŵŝŶĂĐĐĞůĞƌĂƚŝŽŶ͗͞ŽŵĞŐĂͲƐƋƵĂƌĞŵŽĚĞů͘͟ dŚĞƐĞƐƉĞĐƚƌĂĂƉƉĞĂƌ͞ŽŵĞŐĂͲƐƋƵĂƌĞ͟ĨƌŽŵϬ͘ϮͲϭ͘Ϭ,ǌ͘^ŝƚĞZĞƐƉŽŶƐĞ ^ŝŵƉůĞŵŽĚĞů͗ůĂǇĞƌĞĚƐƚƌƵĐƚƵƌĞ͕ŝŶĐůƵĚŝŶŐĂƚƚĞŶƵĂƚŝŽŶ͘ DŽĚĞůĨƌŽŵ<ĂǁĂƐĞ >Ăď͗ŝŶŶŽǀĂƚŝǀĞƵƐĞŽĨŶŽŝƐĞΘĐŽĚĂ͘ ^ƵŐŐĞƐƚƐƚŚĂƚĚĞǀŝĂƚŝŽŶĨƌŽŵ͞ŽŵĞŐĂͲƐƋƵĂƌĞ͟ĂďŽǀĞϭ,ǌĐŽƵůĚďĞĞŶƚŝƌĞůǇĚƵĞƚŽƐŝƚĞƌĞƐƉŽŶƐĞ͘ŽŵƉĂƌĞƐŚĂƉĞǁŝƚŚĚĂƚĂ

>ŽǁĞƐƚĨůĂƌŐĞůǇƵŶĚĞƌƐƚŽŽĚ͕ďƵƚŶŽƚƉƌĞĚŝĐƚĂďůĞĨŽƌĨƵƚƵƌĞƋƵĂŬĞƐ͘>ŽŶŐǁĂǀĞƐʹ ĚŽǁŶŚŽůĞŝƐĂƚ͞ĨƌĞĞƐƵƌĨĂĐĞ͟<ĂǁĂƐĞ ŵŽĚĞů

ŵƉůŝĨŝĐĂƚŝŽŶĚŽŵŝŶĂƚĞƐ

ϯϯ

ϯϰ

−45−

tŚǇĂƌĞŵĞĂƐƵƌĞŵĞŶƚƐŽĨ͞ŬĂƉƉĂ͟ƐŽǀĂƌŝĂďůĞ͍ĂŶƐǇŶƚŚĞƚŝĐƐĞŝƐŵŽŐƌĂŵƐďĞďĞƚƚĞƌƚŚĂŶ'DWƐ͍ <ĂƉƉĂʹ ĂƉĂƌĂŵĞƚĞƌƚŚĂƚĐŚĂƌĂĐƚĞƌŝǌĞƐĂĨŝƌƐƚͲŽƌĚĞƌĨŝƚƚŽƚŚĞƐŚĂƉĞŽĨƚŚĞĂĐĐĞůĞƌĂƚŝŽŶƐƉĞĐƚƌƵŵĂƚŚŝŐŚĨƌĞƋƵĞŶĐŝĞƐ͘ ^ŝŵƉůĞŵŽĚĞůĐĂŶĞǆƉůĂŝŶŝƚƐďĞŚĂǀŝŽƌƵƐŝŶŐǁĂǀĞƉƌŽƉĂŐĂƚŝŽŶĂƚƚĞŶƵĂƚŝŽŶ͘ ,ŽǁĞǀĞƌ͕ƐŽŵĞŽďƐĞƌǀĞĚǀĂƌŝĂďŝůŝƚǇŝƐďĞƚƚĞƌĞǆƉůĂŝŶĞĚƵƐŝŶŐĂƐŽƵƌĐĞƚĞƌŵ͕ƐƵŐŐĞƐƚŝŶŐƚŚĂƚƚŚĞƐŽƵƌĐĞŝƐŶŽƚĐŽŶƐŝƐƚĞŶƚůǇŽŵĞŐĂͲƐƋƵĂƌĞ͘ DĂŶǇĂƉƉƌŽĂĐŚĞƐƚŽƐǇŶƚŚĞƚŝĐƐĞŝƐŵŽŐƌĂŵƐ͘ϯϱ ^ƚŽĐŚĂƐƚŝĐ
 /ƌŝŬƵƌĂ ZĞĐŝƉĞ ŽŵƉŽƐŝƚĞ^ŽƵƌĐĞDŽĚĞů ŽŵďŝŶĞĨŝŶŝƚĞĚŝĨĨĞƌĞŶĐĞŵŽĚĞůƐ;ůŽǁĨͿΘƐƚŽĐŚĂƐƚŝĐ;ŚŝŐŚĨͿ͘ sĞƌǇŝŵƉŽƌƚĂŶƚWZ/ĂĐŚŝĞǀĞŵĞŶƚƐ͘ /ƌŝŬƵƌĂ ƌĞĐŝƉĞʹ ĞŶŽƌŵŽƵƐŝŵƉĂĐƚ͘ϯϲ,ŽǁĞĨĨĞĐƚŝǀĞůǇĐĂŶƐŝƚĞƌĞƐƉŽŶƐĞďĞƌĞƉƌĞƐĞŶƚĞĚŝŶƐǇŶƚŚĞƚŝĐƐĞŝƐŵŽŐƌĂŵƐ͍

6\QWKHWLF6HLVPRJUDPVXVLQJD&RPSRVLWH6RXUFH0RGHO /ůůƵƐƚƌĂƚĞǁŝƚŚĞǆĂŵƉůĞŽĨǁŽƌŬ/͛ǀĞĚŽŶĞǁŝƚŚWƌŽĨ͘<ĂǁĂƐĞ
ĂŶĚŽƚŚĞƌƐŝŶŚŝƐůĂď͘ dŚĞĂƉƉƌŽĂĐŚŝŶƚƌŽĚƵĐĞƐƐŝƚĞƌĞƐƉŽŶƐĞƚŚƌŽƵŐŚĚĞƚĂŝůĞĚŵŽĚĞůŽĨƐŚĂůůŽǁƐƚƌƵĐƚƵƌĞ͘ϯϳϯϴ

'DWD0RGHO)RXULHU6SHFWUXP

−46−
























'DWD0RGHO36$

I+]





 

DŽĚĞůϰ

 

 

 

7V

DŽĚĞůϰ
ϯϵ

ϰϬ

)6









'ĞŶĞƌĂů/ƐƐƵĞƐ

&DVH
&DVH
&DVH
&DVH









 

 





I+]

36$






 

 

7QV
ϰϭ

ϰϮ

−47−
 ƐƐƵŵƉƚŝŽŶƐ͗ ƚƌƵĞ͕ĚĂƚĂͲŐĞŶĞƌĂƚŝŶŐŚĂǌĂƌĚĐƵƌǀĞĞǆŝƐƚƐĨŽƌĞǀĞƌǇƐŝƚĞ͘ ,ŽǁĞǀĞƌ͕ĞƐƚŝŵĂƚĞƐǀĂƌǇĨƌŽŵŽŶĞŵŽĚĞůƚŽĂŶŽƚŚĞƌ͘>ŽŐŝĐdƌĞĞĨŽƌƉŝƐƚĞŵŝĐhŶĐĞƌƚĂŝŶƚǇ
:

'HSWK

,ĂǌĂƌĚŵĂƉƐƐŚŽƵůĚƐŚŽǁƵŶĐĞƌƚĂŝŶƚŝĞƐ͘,ŽǁƐŚŽƵůĚƚŚŽƐĞƵŶĐĞƌƚĂŝŶƚŝĞƐďĞŵĞĂƐƵƌĞĚ͍

(

 7
UX
FN
HH

, QF
OLQ
H



9
LOOD
J

9L
%H 6X
UJ
OLH PP
LQL
YH LW
D&
3 0
OD] DOO
LW\
D 5
5H H Q
QR R

H












 R

 R


'LVWDQFH1RUPDOWR0RXQW5RVH)DXOWNP

 R


 D͗DƵƚƵĂůůǇǆĐůƵƐŝǀĞΘŽůůĞĐƚŝǀĞůǇǆŚĂƵƐƚŝǀĞ͘ tĞŝŐŚƚƐďĂƐĞĚŽŶũƵĚŐĞŵĞŶƚ͕ŶŽƚĞǆƉĞƌŝŵĞŶƚ͘ DƵůƚŝƉůǇǁĞŝŐŚƚƐĨƌŽŵĞĂĐŚŶŽĚĞƚŽŐĞƚďƌĂŶĐŚǁĞŝŐŚƚĂƚƚŚĞĞŶĚ͘ ĂĐŚďƌĂŶĐŚŚĂƐĂĐŽƌƌĞƐƉŽŶĚŝŶŐŚĂǌĂƌĚĐƵƌǀĞ͘ƌĞƚŚĞƐĞĚŝĨĨĞƌĞŶĐĞƐƐŝŐŶŝĨŝĐĂŶƚ͍ϰϯϰϰs^ϯϬ ŝƐĂƉƌŽǆǇ
/ŶE'tĞƐƚϮ&ůĂƚĨŝůĞ

,ŽǁƐŚŽƵůĚƐŝƚĞƌĞƐƉŽŶƐĞďĞƌĞƉƌĞƐĞŶƚĞĚŝŶ
'DWƐ͍ ϯϬǇĞĂƌƐĂŐŽǁĞǁĞƌĞƵƐŝŶŐƋƵĂůŝƚĂƚŝǀĞƚĞƌŵƐƚŽĚĞƐĐƌŝďĞƚŚĞƐŝƚĞƉƌŽƉĞƌƚŝĞƐ͘ ǆƚĞŶƐŝǀĞŵĞĂƐƵƌĞŵĞŶƚƐƐŝŶĐĞƚŚĞŶ͘ ƵƌƌĞŶƚůǇ͕s^ϯϬ ŝƐĂĐŽŵŵŽŶƉƌŽǆǇĨŽƌƐŝƚĞƌĞƐƉŽŶƐĞŝŶ
'DWƐ͘ ĂŶǁĞĚŽďĞƚƚĞƌ͍,ĂƌĚƐŝƚĞƐŽŶƚŽƉ^ŽĨƚƐŝƚĞƐŽŶƚŽƉ

ϮϮ

Ϯϭ

−48−
$SULO)XNXVKLPD+DPDGRULHDUWKTXDNH0:

$SULO)XNXVKLPD+DPDGRUL HDUWKTXDNH0:
,PDJHVIURP0DUX\DPDHWDO *6-

ĂďŝŶĂĚũĂĐĞŶƚƚŽĨĂƵůƚƚŚĂƚƌƵƉƚƵƌĞĚŝŶƚŚĞϭϵϱϰŝǆŝĞsĂůůĞǇĞĂƌƚŚƋƵĂŬĞ;Dtϲ͘ϵͿ͘dŚĞƉŚŽƚŽŐƌĂƉŚǁĂƐƚĂŬĞŶŽŶϮϬĞĐĞŵďĞƌϭϵϱϰďǇ<Ăƌů^ƚĞŝŶďƌƵŐŐĞ ;E/^͕hͿ͘0DUFK0:,PPHGLDWHDIWHUVKRFNV&ĂƵůƚĞǆƉŽƐƵƌĞĚƵƌŝŶŐƚŚĞϮϬϭϭ&ƵŬƵƐŚŝŵĂ,ĂŵĂĚŽƌŝ ĞĂƌƚŚƋƵĂŬĞ͘dŚĞƐĞƉŚŽƚŽŐƌĂƉŚƐǁĞƌĞƚĂŬĞŶďǇdĂĚĂƐŚŝDĂƌƵǇĂŵĂ͕'ĞŽůŽŐŝĐĂů^ƵƌǀĞǇŽĨ:ĂƉĂŶ͘ϮϯϮϰ^ƵŵŵĂƌǇ͘DǇůŝƐƚŽĨŬĞǇƋƵĞƐƚŝŽŶƐĨŽƌƐĞŝƐŵŝĐŚĂǌĂƌĚĞƐƚŝŵĂƚŝŽŶ͘^ƵŵŵĂƌǇ͘DǇůŝƐƚŽĨŬĞǇƋƵĞƐƚŝŽŶƐĨŽƌƐĞŝƐŵŝĐŚĂǌĂƌĚĞƐƚŝŵĂƚŝŽŶ͘ϭ͘ ^ŚŽƵůĚƚŚĞƐĞŝƐŵŝĐŝƚǇĐĂƚĂůŽŐďĞĚĞĐůƵƐƚĞƌĞĚ͍/ĨƐŽ͕ŚŽǁ͍Ϯ͘ tŚĂƚŝƐƚŚĞŵĂŐŶŝƚƵĚĞʹ ĨƌĞƋƵĞŶĐǇĚŝƐƚƌŝďƵƚŝŽŶĨŽƌƚŚĞŶĞƚǁŽƌŬŽĨĨĂƵůƚƐ͍ϯ͘ ,ŽǁĚĞĞƉůǇĚŽůĂƌŐĞĐŽŶƚŝŶĞŶƚĂůĞĂƌƚŚƋƵĂŬĞƐƌƵƉƚƵƌĞ͍/ƐƚŚĞĚĞƉƚŚƉƌĞĚŝĐƚĂďůĞ͍ϰ͘ /ƐŵŽŵĞŶƚďĂůĂŶĐŝŶŐďĂƐĞĚŽŶƐůŝƉͬƐƚƌĂŝŶƌĂƚĞƚŚĞďĞƐƚĂƉƉƌŽĂĐŚƚŽĞƐƚŝŵĂƚĞĞĂƌƚŚƋƵĂŬĞƌĂƚĞƐ͍ƌĞƚŚĞƌĞŽƚŚĞƌŵĞƚŚŽĚƐ͍ϱ͘ tŚĂƚŝƐƚŚĞŐƌŽƵŶĚŵŽƚŝŽŶŶĞĂƌĂŐƌĞĂƚĞĂƌƚŚƋƵĂŬĞŝŶƚŚĞĐŽŶƚŝŶĞŶƚĂůĐƌƵƐƚ͍ϱ͘ tŚĂƚŝƐƚŚĞŐƌŽƵŶĚŵŽƚŝŽŶŶĞĂƌĂŐƌĞĂƚĞĂƌƚŚƋƵĂŬĞŝŶƚŚĞĐŽŶƚŝŶĞŶƚĂůĐƌƵƐƚ͍
ϲ͘ ,ŽǁƐŚŽƵůĚƐŝƚĞƌĞƐƉŽŶƐĞďǇƌĞƉƌĞƐĞŶƚĞĚŝŶ'DWƐ͍ϳ͘ ĂŶƐŝŐŵĂďĞƌĞĚƵĐĞĚŝŶW^,͍ϴ͘ ,ŽǁĚŽǁĞĂĚũƵƐƚŐƌŽƵŶĚŵŽƚŝŽŶŵŽĚĞůƐĨŽƌĚŝĨĨĞƌĞŶƚƌĞŐŝŽŶĂůŐĞŽůŽŐǇĂŶĚƚĞĐƚŽŶŝĐƐ͍ϵ͘ tŚĂƚŝƐƚŚĞƐŚĂƉĞŽĨƚŚĞƐĞŝƐŵŝĐ͞ƐŽƵƌĐĞƐƉĞĐƚƌƵŵ͕͟;ĨͿ͍ϭϬ͘ tŚǇĂƌĞŵĞĂƐƵƌĞŵĞŶƚƐŽĨ͞ŬĂƉƉĂ͟ƐŽǀĂƌŝĂďůĞ͍ϭϭ͘ ĂŶƐǇŶƚŚĞƚŝĐƐĞŝƐŵŽŐƌĂŵƐďĞďĞƚƚĞƌƚŚĂŶ'DWƐ͍ϭϮ͘ ,ŽǁĞĨĨĞĐƚŝǀĞůǇĐĂŶƐŝƚĞƌĞƐƉŽŶƐĞďǇƌĞƉƌĞƐĞŶƚĞĚŝŶƐǇŶƚŚĞƚŝĐƐĞŝƐŵŽŐƌĂŵƐ͍

ϰϵ

ϱϬ

−49−

^ƵŵŵĂƌǇ͘DǇůŝƐƚŽĨŬĞǇƋƵĞƐƚŝŽŶƐĨŽƌƐĞŝƐŵŝĐŚĂǌĂƌĚĞƐƚŝŵĂƚŝŽŶ͘tĞŶĞĞĚƚŽǁŽƌŬƚŽŐĞƚŚĞƌ͊ϭϯ͘ ,ĂǌĂƌĚŵĂƉƐƐŚŽƵůĚƐŚŽǁƵŶĐĞƌƚĂŝŶƚŝĞƐ͘,ŽǁƐŚŽƵůĚƚŚŽƐĞƵŶĐĞƌƚĂŝŶƚŝĞƐďĞŵĞĂƐƵƌĞĚ͍ϭϰ͘ ,ŽǁĐĂŶŚĂǌĂƌĚŵĂƉƐĂŶĚŚĂǌĂƌĚĐƵƌǀĞƐďĞƚĞƐƚĞĚ͍ϭϱ͘ ŽǁĞŶĞĞĚƚŽĐŚĂŶŐĞŽƵƌƉƌĂĐƚŝĐĞƐƚŽĂĐŚŝĞǀĞŐƌĞĂƚĞƌƌĞƐŝůŝĞŶĐĞ͍ dŚĂŶŬǇŽƵ͘ &ŽƌǇŽƵƌĂƚƚĞŶƚŝŽŶƚŽĚĂǇ
 &ŽƌƚŚŝƐŐĞŶĞƌŽƵƐƌĞĐŽŐŶŝƚŝŽŶ͘ &ŽƌǁĞůĐŽŵŝŶŐŵĞƚŽWZ/ƐŽŵĂŶǇƚŝŵĞƐŝŶƚŚĞƉĂƐƚ͘ &ŽƌĐŽůůĂďŽƌĂƚŝŽŶĂŶĚĨƌŝĞŶĚƐŚŝƉ͘ϱϭϱϮ−51−

−52−

−53−

−54−

−55−

−56−

−57−

−58−

−59−

−60−

−61−

−62−

−63−

−64−

−65−

−66−

−67−

−68−

−69−

−70−

−71−

−72−

−73−

−74−

−75−

−77−'RZQORDGHGIURPZZZVFLHQFHPDJRUJRQ2FWREHU−78−'RZQORDGHGIURPZZZVFLHQFHPDJRUJRQ2FWREHU−79−'RZQORDGHGIURPZZZVFLHQFHPDJRUJRQ2FWREHU−80−'RZQORDGHGIURPZZZVFLHQFHPDJRUJRQ2FWREHU−81−'RZQORDGHGIURPZZZVFLHQFHPDJRUJRQ2FWREHU−82−'RZQORDGHGIURPZZZVFLHQFHPDJRUJRQ2FWREHU−83−'RZQORDGHGIURPZZZVFLHQFHPDJRUJRQ2FWREHU

Bulletin of the Seismological Society of America, Vol. 103, No. 3, pp. 1935–1951, June 2013, doi: 10.1785/0120120252

Ground Motions in the Fukushima Hamadori,
Japan, Normal-Faulting Earthquake
by John G. Anderson,* Hiroshi Kawase, Glenn P. Biasi, James N. Brune, and Shin Aoi
A crustal normal-faulting earthquake (M JMA 7:0; M w 6.7) occurred ineastern Tohoku, Japan, on 11 April 2011. K-NET and KiK-net stations recorded
82 records from within 100 km of fault rupture. These data and data from associatedforeshocks and aftershocks will make a critical contribution to future improvements ofground-motion prediction for normal-faulting earthquakes.
Peak ground accelerations (PGA) and peak ground velocities (PGV) are comparedwith four ground-motion prediction equations (GMPEs) that include the style of faulting as a predictor parameter. For distances under 100 km, and using a network averagevalue of V S30 , the average ratio of PGA to the selected GMPEs (the event term ) is highby factors of 2.3–3.7. Event terms for PGV are high by factors of 1.4–1.8. Adjusting
PGA and PGV with customized site terms (Kawase and Matsuo, 2004a,b), the standard
deviations of PGA and PGV residuals are reduced from 0.59 to 0.43, and from 0.53 to
0.35, respectively. The event terms decreased to relatively small factors of 1.1–1.8 for
PGA and increased slightly to 1.5–2.0 for PGV. Thus, site terms are very important, but
positive event terms remain. The remaining positive event terms are not explained byhigh stress drop, which was typical of crustal events of all mechanisms globally or inJapan. Two subparallel faults ruptured, but source inversions, which we reviewed,
revealed that they ruptured sequentially, so simultaneous contributions from the two
faults did not cause high motions. Although these observations may tend to suggestthat ground motions in large normal-faulting events are larger than predicted by thetested models, we are not aware of any observations from this event that contradict the
precarious rock evidence of Brune (2000) that ground shaking is low on the footwall
near the rupture.

Abstract

Introduction
Strong-motion records from the near field of earthquakes with a normal-faulting mechanism (Fig. 1) are sparse,
even by comparison with crustal earthquakes with a strikeslip or reverse-faulting mechanism. For instance, a search of
the Next Generation Attenuation (NGA) relation flatfile(Chiou et al., 2008) finds that only 91 records (out of 3551
total used for the NGA) have a normal-faulting mechanism.
There are only 25 records from three normal-faulting earthquakes with M w ≥ 6:5 (Irpinia, Italy, 23 November 1980,
Mw 6.9, 12 records; Borah Peak, Idaho, 20 October 1983,
Mw 6.9, 11 records [in a fairly small area]; Edgecumbe, New
Zealand, 2 March 1987, M w 6:6, 2 records). The NGA compilation, of course, does not include the entire set of groundmotions from normal-faulting earthquakes. It has the obvious
limitation that the most recent event in that compilation occurred in 2003. Furthermore, it did not benefit from recent*Also at Department of Geological Sciences and Engineering.

efforts to improve the quality and accessibility of strongmotion data from Italy (e.g., Luzi et al., 2010; Massa et al.,2010; Pacor et al., 2011) and Greece (Theodulidis et al.,2004), where normal faulting is dominant.
Within these modest constraints, several ground-motion
prediction equations (GMPEs) have concluded that on average, ground motions in normal-faulting events are slightly
smaller than in strike-slip or thrust events (e.g., Spudich et al.,
1999; Abrahamson and Silva, 2008; Boore and Atkinson,
2008; Campbell and Bozorgnia, 2008; Chiou and Youngs,
2008). The last four of these are collectively known as NGA
relations. A recent model by Bindi et al. (2011) using the
Italian database finds that ground motions in Italy are generally consistent with the NGA predictions for normal-faulting
events, but near-source records from earthquakes with
Mw ≥ 6 are sparse in that data set as well. The NGA relations,developed under the NGA Project (Power et al., 2008), areused extensively in hazard analysis in the United States.
Applications include use of Boore and Atkinson (2008),

1935

−85−

1936
7

J. G. Anderson, H. Kawase, G. P. Biasi, J. N. Brune, and S. Aoi
Normal Faulting Earthquakes in the NGA Flatfile

731

40.0° N

Magnitude

6.5

6
516
393

5.5

454
415
350

37.5° N

320
244
258

5

4.5
0
10

207

1

2

10

10

3

10

r , kmflt

35.0° N

Figure 1.

Distribution of magnitude and distance of strongmotion data for normal faulting earthquakes in the NGA database,as a function of rflt . For events in which rflt is not given in the database, repi is substituted. The solid horizontal lines are drawn at magnitudes of normal faulting earthquakes associated with the
Fukushima Hamadori earthquake. The numbers to the right givethe number of strong-motion accelerograms recorded in each event,and the extent of the lines give a preliminary estimate of the distance
range of each event. Distances are based on epicenter as given in the
K-NET or KiK-net files obtained from NIED.

Campbell and Bozorgnia (2008), and Chiou and Youngs(2008) by the U.S. Geological Survey to develop the 2008
revision of the National Hazard Map (Petersen et al., 2008).
These three models are used to estimate the ground motions
from earthquakes in the Intermountain West region, where
normal faulting dominates.
In this context, the Fukushima Hamadori earthquake of
11 April 2011 (M w 6.7, MJMA 7.0) is extremely significant,as the strong ground motions are by far the most extensively
recorded of any normal-faulting earthquake of any magnitude. About 82 stations are within 100 km of the nearer of
the two fault surfaces, at a wide range of azimuths, and over
730 stations altogether recorded the event. Besides the mainshock, the Global Centroid Moment Tensor (Global CMT)
project (see Data and Resources) identified at least nine othernormal-faulting events in this sequence with M w > 5, as indicated in Figure 1. This analysis, which evaluates only themainshock, finds that the recorded strong motion in the
Fukushima Hamadori mainshock had significantly stronger
peak ground accelerations (PGA) and peak ground velocities(PGV) than the prediction of the NGA equations cited above,and considers four hypotheses that might help explain thediscrepancy.

Fukushima Hamadori Earthquake
The Fukushima Hamadori earthquake, and associatedevents with smaller magnitudes, was apparently triggeredby the 11 March 2011 Tohoku (M w 9.0) earthquake. The
Fukushima Hamadori earthquake occurred in the continental

140.0° E

142.5° E

145.0° E

Figure 2.

Location map, showing the region affected by the
Fukushima Hamadori earthquake (outlined by heavy box) in thecontext of the region affected by the Tohoku earthquake. Epicentersshown are events located by the U.S. Geological Survey National
Earthquake Information Center (USGS NEIC) in the first 12 days
following the M w 9.0 Tohoku mainshock. Dashed black lines indicate the boundary between the Pacific or Philippine Sea plates andthe crustal block including Tohoku.

crust of Japan, and caused surface rupture. Large shallow
aftershocks with extensional mechanisms also followed theMw 8.8 Maule, Chile, earthquake in 2010 (Ryder et al., 2012).
The general location of the Fukushima Hamadori event, in thecontext of the Tohoku aftershocks, is shown in Figure 2, andsome teleseismic parameters for the earthquake are given in
Table 1. Although Kato et al. (2011) suggested that the entire
Tohoku region, including the source area of the Fukushima
Hamadori earthquake, was under compression prior to the
Tohoku earthquake, Imanishi et al. (2012) found 23 focal
mechanisms of smaller earthquakes with normal-faulting
mechanisms in the region between 2003 and 2010.
The Fukushima Hamadori earthquake caused surface
rupture on two faults. Faulting is described by Maruyama
et al. (2011). The Idosawa fault ruptured for a distance of
14 km, with a maximum offset of about 200 cm (Maruyama
et al., 2011). The Yunotake fault ruptured for about 15 km,
with maximum slip of about 80 cm. Table 2 gives the parameters used in this paper to describe the locations and orientations of these two faults. A model of the faulting is needed
to determine the fault distance for comparisons of observed
ground motions with GMPEs. Figure 3 shows the locationof the fault planes as modeled in Table 2. Both faults areapproximated by a planar surface, with the surface trace
approximating the location of observed rupture as given by
Maruyama et al. (2011) and Geology Research Team (2011).
The dip of the Idosawa fault is constrained so that the fault−86−Ground Motions in the Fukushima Hamadori, Japan, Normal-Faulting Earthquake

1937

Table 1
Catalog Parameters for the Fukushima Hamadori Earthquake of 11 April 2011
Data Source

Origin Time

K-NET and KiK-net
data files

Depth, km

M 0 dyn·cm

Magnitude

6.0

7:0M JMA

08:16:12.73 UTC

11.0

08:16:19.5 UTC (centroid time)

37.00° N
140.54° E

12.0

6:6Mw W
6:69Mw
6:4mb
6:6MS

17:16 JST

PDE

Global CMT *Location

36.945° N
140.672° E
37.001° N
140.401° E

1:2 × 1026

*Fault planes: FP1: Strike  139, Dip  56, Rake  −75; FP2: Strike  295, Dip  39, Rake  −109; FP1 ismost consistent with the observed surface faulting.

Strong-Motion Data

Table 2
Geometry of Fault Planes Used in This Study to
Describe the Idosawa and Yunotake Faults
Parameter

Idosawa

Yunotake

Reference point latitude*
Reference point longitude*
Reference point depth (km)*Strike†Dip†
Total Rupture Length (km)†Total Rupture Width (km)†Length in +strike direction fromreference point (km)‡Length in −strike direction fromreference point (km)‡Shallowest depth of faulting (km)Maximum depth of faulting (km)
Fault area (km2 )Maximum Surface Slip (m)36.9488
140.7082
0
162
66
13.4
20
0

36.9855
140.8302
0
128
66
15.3
15
0

13.4

15.3

0
18.27
268
2.0

0
13.70
229.5
0.8

*The reference points are at the south end of the surface traces ofeach fault.
†The strike and fault length are chosen to be a close approximation to the location of surface faulting. The strike uses the Aki
and Richards (2002) convention that when facing in the +strikedirection, the fault dips to the right. The dip of the Idosawa fault
is determined so that the dipping plane will go through the hypocenter, whereas the dip of the Yunotake fault is taken to be the same.
The down-dip width is wide enough to reach the depth and horizontal
offset of the deepest aftershocks. Cross sections of the faults showthat these planes are not inconsistent with the aftershock locations.
‡Because of the convention for fault geometry and the locationof the reference points, the fault traces are in the opposite direction
from the positive-strike direction.

will pass through the hypocenter, and the same dip is thenapplied to the Yunotake fault. Associated smaller events thatwere also recorded on strong-motion instruments include immediate aftershocks, and the fault width in Table 2 was set
so as to be consistent with the deepest of those events. Basedon aftershocks, Hikima (2012) and Shiba and Noguchi(2012) model the Idosawa fault as longer than its surface rupture, 26 and 23.8 km, respectively, mainly extending the fault
towards the north compared with Figure 3. The kinematic
models by both authors find most of the slip took place
within the part of the fault shown in the figure.

Strong-motion data in this study was obtained from the
NIED K-NET and KiK-net strong-motion networks (Kinoshita, 1998, 2005; Okada et al., 2004). Data from all stations
recording the event were downloaded from the NIED websites. There are 82 stations within 100 km of the faults shownin Figure 3, 298 within 200 km, 477 stations within 300 km,
608 stations within 400 km, and 715 stations altogether, withthe most distant at 876 km distance. Data were high-passfiltered with a Butterworth causal filter, 2 pole, with cornerfrequency of 0.05 Hz, and then integrated to velocity anddisplacement. The geometric mean horizontal accelerationand geometric mean horizontal velocity, determined from thehorizontal components oriented as recorded, are shown as afunction of distance in Figures 4 and 5, respectively. FourNGA models (Abrahamson and Silva, 2008; Boore and Atkinson, 2008; Campbell and Bozorgnia, 2008; Chiou andYoungs, 2008) are also shown in Figures 4 and 5. These four
NGA relations use V S30 , the time-average shear velocity inthe upper 30 m, as a predictor variable. The stations contributing data on Figures 4 and 5, of course, come from a variety
of site conditions. Therefore, the intention on these figures is
to use a value of V S30 that is a reasonable median for theensemble of stations. All of the K-NET sites have velocity
logs posted by NIED, giving compressional and shear velocities and densities at one meter intervals from the surface to
20 m depth. A reasonable lower bound to V S30 , say V min
S30, is toproject the velocity at the base of the provided log to thedepth of 30 m. The median value of V min
S30 is 340 m=s. We also,
byprojecting a velocity
estimated an upper bound, say V max
S30
of 2000 m=s from the base of velocity log to 30 m. The
median value of V max
S30 estimated in this manner is 470 m=s.
Therefore, the median V S30 is bounded between 340 and
470 m=s. The value V S30  400 m=s is used in Figures 4
and 5. It is worth noting that based on the correlation of
Boore et al. (2011), a median value of V S20  340 m=s corresponds to a median value of V S30  407 m=s.
The data are offset above the predictions in Figures 4
and 5. A different choice of V S30 can not help significantly
to reduce the discrepancy. For several spot checks forboth PGA and PGV at distances of 10, 50, and 100 km,

−87−

1938

J. G. Anderson, H. Kawase, G. P. Biasi, J. N. Brune, and S. Aoi
3

10

PGA (Mean Horizontal), cm/s2

37.5° N

37.0° N

2

10

1

10

0

10

−1

Yunotake
fault

Data
Abrahamson and Silva (2008)
Boore and Atkinson (2008)
Campbell and Bozorgnia (2008)
Chiou and Youngs (2008)

10

1

2

10

10

3

10

rflt, km

Idosawa fault

Figure 4. PGA, compared with four different GMPEs. Curves
assume V S30  400 m=s, as explained in the text.

2

10

140.0° E

140.5° E

141.0° E

Figure 3.

Detailed map of the area around the Fukushima Hamadori earthquake. Black rectangles show the surface projection of
the two contributing faults, as modeled in Table 2. The fault edges
with tick marks indicate the approximate locations of surface rupture. Triangles show locations of strong-motion stations. The starshows the epicenter of the mainshock. Circles are associated crustal
earthquakes, which have generated strong-motion accelerograms.
Earthquake locations are NIED locations from the headers ofstrong-motion records obtained at station IBRH13 from shallow,
nearby earthquakes.

the difference between curves generated with V S30 
400 m=s and the predictions at either bound on median V S30
is under 6%. This difference is negligible compared with theoffset of the data.
The GMPE model by Si and Midorikawa (1999, 2000) is
used in the development of the national seismic hazard mapof Japan. It is not shown in Figures 4 and 5 because thatmodel is developed for PGA and PGV on records that havebeen low-pass filtered with a corner frequency of 10 Hz,
and it predicts the value of the larger horizontal component.
In contrast, the accelerograms used to develop Figures 4 and5 show the geometric mean of the two horizontal peaks without applying a low-pass filter to the data. The Si and Midorikawa (1999, 2000) model predicts larger amplitudes of
PGA and PGV, for distances smaller than 100 km, than the
NGA models, but the observations of both PGA and PGV are
still somewhat higher (Si et al., 2012). Because the Si andMidorikawa (1999, 2000) model does not include an adjustment for fault mechanism, this implies that ground motionsin the Fukushima Hamadori earthquake are somewhat
stronger than the average reverse and strike-slip earthquakes

PGV (Mean Horizontal), cm/s

36.5° N

1

10

0

10

−1

10

−2

10

Data
Abrahamson and Silva (2008)
Boore and Atkinson (2008)
Campbell and Bozorgnia (2008)
Chiou and Youngs (2008)
1

10

2

10

3

10

rflt, kmFigure 5.

PGV compared with four different GMPEs. Curves assume V S30  400 m=s, as explained in the text.

that dominate in development of the Si and Midorikawa(1999, 2000) model.

Event Terms
The value of an individual ground-motion observation,using the notation of Anderson and Uchiyama (2011), can be
written asln Y es  fM ·  δ Ies ;(1)

where Y es is the observation (either PGA or PGV), fM · isthe value predicted by the GMPE, and δ Ies is the residual. The
subscripts e and s identify the earthquake (just one in thiscase) and the station. The function fM · will be different for−88−

Ground Motions in the Fukushima Hamadori, Japan, Normal-Faulting Earthquake

1939

Table 3
Event Terms for the NGA GMPEs Determined from Stations at Distances under 100 km
Acceleration

Velocity

GMPE

Ee

εa

σIIE

e Ee

Ee

εv

σIIE

e Ee

Abrahamson and Silva (2008)
σIa  0:573, σIv  0:583
σEvents;a  0:326, σEvents;v  0:319
Boore and Atkinson (2008)
σIa  0:564, σIv  0:560
σEvents;a  0:260, σEvents;v  0:256
Campbell and Bozorgnia (2008)
σIa  0:515, σIv  0:525
σEvents;a  0:219, σEvents;v  0:203
Chiou and Youngs (2008)
σIa  0:549, σIv  0:540
σEvents;a  0:271, σEvents;v  0:241
Average

1.08

3.3

0.59

2.95

0.60

1.9

0.52

1.82

0.84

3.2

0.57

2.32

0.51

2.0

0.53

1.67

0.98

4.5

0.59

2.66

0.32

1.6

0.52

1.37

1.32

4.8

0.62

3.74

0.54

2.2

0.54

1.71

2.87

0.49

1.05

PGA and PGV. Higher roman numeral superscripts are used
to indicate standard deviations after removal of potential
epistemic terms. The standard deviation corresponding toequation (1) iss
1 X I 2δ  ;σ 
N es − 1 e;s esI(2)

where the total number of data is N es . The event term Ee is
the mean of the residuals in a single event over all stations:
s
1 X
δI ;
N s s1 es

N

Ee 

(3)

where N s is the number of stations observing the event. Asa measure of the variability of events used to develop the
GMPE, it is convenient to define σEvents as the standard
deviation of Ee , and the deviation of a single event from average as
Ee  εσEvents ;(4)(in the notation of Rodriguez-Marek et al., 2011, the event
term Ee is the between-event residual or interevent residual,and σEvents  τ ). Because the main hazard is at close distances, Table 3 gives Ee for PGA and PGV, for stations at
distances under 100 km. GMPEs considered are the NGA relations of Abrahamson and Silva (2008), Chiou and Youngs(2008), and Campbell and Bozorgnia (2008) at rflt < 100 km,
and for Boore and Atkinson (2008) with rJB < 100 km. The
I
unexplained residual, δ IIE
es  δ es − Ee (δW es in the notation ofRodriguez-Marek et al., 2011) has zero mean and standard
deviation σIIE (ϕ in the notation of Rodriguez-Marek et al.,2011), also given in Table 3. The σIIE corresponds to thewithin-event standard deviations of the GMPEs. These σIIE
are not optimized because the site conditions are not treatedon a station-by-station basis.

1.63

The GMPE column in Table 3 gives the average standard
deviation, σI for PGA and PGV, and σEvents , the standard
deviation of the event terms as determined by each respectivemodel. For models in which the standard deviation depends
on distance and amplitude, σI is the average of the value appropriate for each datum. This is a reasonable approximation
because the distance dependence of σI is weak, if present, inthe models. The columns with εa and εv in Table 3 are theevent terms normalized by σEvents as in equation 4. For instance, the event term for PGA relative to the Abrahamson
and Silva (2008) relationship is 3:3 standard deviationsof the event term, determined from all earthquakes used
to develop the model, above the median curve. The columns
with eE convert the log values of the event terms into numerical ratios. Thus, the mean PGA observation is 2.95 times theprediction of Abrahamson and Silva (2008).
Table 3 finds that for PGA the range of σIIE , from 0.57–
0.62, is similar to that found by Anderson and Uchiyama(2011) in Guerrero (0.61), similar to the estimate of Rodriguez-Marek et al. (2011) for Japan (0.62), and similar to thestandard deviations of the NGA models. For PGV the range ofσIIE from 0.52–0.54 is slightly smaller than the Guerrero estimate (0.57), but similar to the standard deviations of the
NGA models. Thus, the variability of these uncorrected data
are in a reasonable range.
On the other hand, Table 3 finds that the event terms arehigh for all four NGA models. PGAs are 3.2–4.8 standard deviations above the NGA models, with an average of about
4 times greater than the respective models (∼  4σEvent;a ).
PGVs average about 2 times greater (∼  2σEvent ). The deviations of PGV would have a low probability, but the deviations for PGA would seem to have such a low probability as
to demand further study.

Discussion: What Causes the High Event Terms?
We consider four possible explanations for the highevent terms of this earthquake at short distances: (1) site

−89−

1940

J. G. Anderson, H. Kawase, G. P. Biasi, J. N. Brune, and S. Aoi

Figure 6.

KiK-net maps showing (a) station name, (b) peak surface/downhole spectral ratio, (c) surface PGA (cm=s2 ), and (d) surface
PGV (cm=s). Peaks of acceleration and velocity are the peak amplitude of the vector sum of the two horizontal components. The color versionof this figure is available only in the electronic edition.

effects are inadequately explained using V S30 ; (2) the source
has a high stress drop, or is otherwise radiating excessive
energy at the frequencies controlling PGA; (3) the two subparallel faults are radiating simultaneously, increasing theamplitudes of shaking; or (4) the GMPEs are underpredicting
the ground motions. Directivity was not considered because
the source is surrounded by stations. The different slopes of
the NGA curves and the data beyond 100 km (Figs. 4 and 5),
probably result from differences in attenuation but large
event terms at short distances are difficult to explain by differences in attenuation.
Site Effects
The comparisons in Figures 4 and 5 already incorporateparameter choices in the NGA relations appropriate for the
median value of V S30 . However, there are reasons to expectthat V S30 by itself is not sufficient to characterize the site
conditions. For instance, Atkinson and Casey (2003) compared amplitudes of pseudoacceleration response spectra(5% damping) at four oscillator frequencies from earth-quakes in southern Japan and Washington State of similar
mechanism and magnitude. Noting that different velocity
profiles can have the same value of V S30 , they explained
higher average amplitudes at 10 Hz in the Japanese eventby the differences in average site profiles. Amplification
was calculated using the quarter-wavelength method (Boore
and Joyner, 1997). This method does not recognize resonantpeaks, but for constant V S30 a site with a thin low-velocity
layer will have higher amplification than a site with a lowvelocity gradient. Thus, characterizing the site response using
V S30 alone is inadequate to fully characterize site response.
Before adjusting for site effects, we view some seismograms and their spectra. Figure 6 shows maps of KiK-net
stations within about 50 km of the Fukushima Hamadori
source. All stations in KiK-net have downhole accelerometers at a typical depth of 100 m in addition to surface
accelerometers. Although the downhole instrument may beaffected by site response to some extent, the differences
between the surface and downhole records provide much information about site response at these stations. Figure 6 also
maps the peak vector acceleration, peak vector velocity, and−90−

Ground Motions in the Fukushima Hamadori, Japan, Normal-Faulting Earthquake

1941

Figure 7.

Velocity seismograms and Fourier amplitude spectra at two KiK-net stations. Station locations are shown in Figure 6. The
velocity seismograms are from the surface components. For FKSH14, rflt  13:4 km, aGM  129 cm=s2 , and vGM  23:1 cm=s. For
IBRH12, rflt  33:5 km, aGM  208 cm=s2 , and vGM  6:3 cm=s. Fourier spectra show Fourier amplitudes of acceleration. Spectra of
the two horizontal components are smoothed in a way that approximately conserves the energy (or to be more precise, the Arias intensity)of the spectrum, and then the vector sum of the two horizontal components is formed. Spectra processed in this way are shown for both thesurface and downhole components, along with the ratio of surface/downhole.

the peak of the ratio of smoothed surface to downhole spectral amplitudes. A peak spectral ratio of under 6.0 was usedby Anderson (2013) to recognize sites with a minimum distortion from surface layers. Figure 6 indicates that none of
the sites in this region meet that criterion for the Fukushima
Hamadori earthquake, but some are close.
Figure 7 shows velocity seismograms and Fourier spectra of acceleration for two of the KiK-net sites. Station
FKSH14 has one of the lower surface/downhole ratios on
Figure 6b. The surface spectrum is relatively flat to first order, with some fine structure superimposed on the overall
spectral shape. In contrast, station IBRH12 has a prominent
narrow peak at about 5–7 Hz. The high peak is associatedwith conspicuous ringing that dominates the surface accelerogram and is even conspicuous on the surface-velocity
seismograms. This record contrasts sharply with the record

from station FKSH14 in which the velocity seismogramshows relatively very little high-frequency shaking. Inspection of seismograms and spectra at the KiK-net sites suggestthat IBRH12 is more typical of accelerograms from thisregion. Stations with high surface/downhole spectral ratios
in Figure 6b show similar phenomena to the resonance at
IBRH12, and contribute high positive residuals in Figure 4.
Figure 8 shows maps of K-NET stations, at which there
are only surface accelerometers. Instead of spectral ratios,
this figure shows V min
S30 for each station. Perusal of Figure 8
shows that there is a correlation between high values of V min
S30
and low PGA. Some of the Fourier spectra show prominent
peaks like the peak at station IBRH12 in Figure 7; others do
not. Figure 9 shows velocity seismograms and Fourier spectra for two stations, FKS011 and FKS014, which do not haveprominent resonances. The residuals at stations FKSH14,

−91−

1942

J. G. Anderson, H. Kawase, G. P. Biasi, J. N. Brune, and S. Aoi

Figure 8.

K-NET composite figure showing (a) station name, (b) estimates of V S30 in m=s, (c) PGA in cm=s2 , and (d) PGV in cm=s. Peaksof acceleration and velocity are the peak amplitude of the vector sum of the two horizontal components. The estimates of V S30 are lower
bounds based on profiles provided by NIED. The site profiles are available to depths of up to 20 m, so these estimates extend the deepest
measured value to 30 m. Any station, for which velocity increases between the bottom of the borehole and 30 m, will have a larger value of
V S30 than the value given here. The color version of this figure is available only in the electronic edition.

FKS011, and FKS014 are all reasonably consistent with the
GMPEs.
First Adjustment for Site Response. The results shown in
Figures 7 and 9 suggest that, as in the case of Atkinson
and Casey (2003), a systematic adjustment for site response
at all the stations might reduce or even eliminate the discrepancy between the data and models. They further suggest thatthe ground motions on the most rock-like stations might berelatively consistent with the NGA relations.
Several models have already been developed to characterize the site response, including nonlinear behavior, at allor most of the individual K-NET and KiK-net stations (e.g.,
Moya and Irikura, 2003; Kawase and Matsuo, 2004a,b;
Wu et al., 2009). In this section, we therefore apply one ofthese models that applies station-specific adjustments foundusing prior earthquakes to reduce each of the KiK-net and
K-NET peak motion observations in the Fukushima Hamadori earthquake to motion expected on rock site conditions.
This differs from Atkinson and Casey (2003), who used corrections for large site classes rather than using site factorsthat are customized on a station-by-station basis.
The model we use is described in Kawase and Matsuo(2004a,b) and Kawase (2006). Kawase and Matsuo (2004a)

defined a reference site with shear velocity of 3:45 km=s that
extends to the surface. Then using moderate-sized earthquakes for each station separately, they found correction factors that adjust the observed PGA and PGV of a record to thepeak acceleration and velocity at the reference site with theshear velocity of 3:45 km=s.
GMPEs used in this paper are defined as a function of
V S30 , but V S30  3450 m=s is outside of their ranges of validity. Thus, we select V S30  760 m=s for the comparison,because that is the reference site condition for the U.S.
National Hazard Map. To adjust the difference of a referencesite between V S30  3450 m=s and V S30  760 m=s, we use
the empirical relationship for average amplification as a loglinear function of V S30 (Kawase and Matsuo, 2004b). Thus,
the observed value of PGA and PGV is divided by the sitespecific adjustment, and then multiplied by the reference correction factor (1.91 for PGA and 2.66 for PGV).
Figure 10 shows adjustments from raw PGA and PGV toa site with V S30  760 m=s, as a function of estimates of
V S30 . The estimates of V S30 are not used in our procedure,
but they provide a convenient parameter to use for displaying
the adjustments. For K-NET stations, Figure 10 uses V min
S30 .
For KiK-net stations, the NIED website gives measured values of shear velocity to depths of 100 m or more that are used−92−

Ground Motions in the Fukushima Hamadori, Japan, Normal-Faulting Earthquake

1943

Figure 9. Velocity seismograms and Fourier amplitude spectra at two K-NET stations. Station locations are shown in Figure 8. For
station FKS011, rflt  13:3 km, aGM  188 cm=s2 , and vGM  17:6 cm=s and for station FKS014, rflt  24:0 km, aGM  111 cm=s2 ,and vGM  8:7 cm=s. Fourier spectra are smoothed as in Figure 7.
to calculate V S30 . Figure 10 also shows adjustments used by
Atkinson and Casey (2003), which do depend on estimatesof V S30 . Notable features of Figure 10 are the high variability of the adjustments within small ranges of V S30 , the tendency of adjustments to decrease as V S30 increases, and thatthe adjustments appear on average to be smaller than those of
Atkinson and Casey (2003).
Although the adjustments in Figure 10 are determinedfrom a different data set, this step can be considered approximately equivalent to the station correction step in Andersonand Uchiyama (2011). Designating the Kawase and Matsuo(2004a,b) corrections as Ss (δS2Ss in the notation of Rodriguez-Marek et al., 2011), the new residuals areIδIIS
es  δ es − Ss :(5)

New estimates of the event terms (E2
e ) can be found by averaging δIIS
es . The residuals after removing these event terms are

2
IIS
δ III
es  δ es − Ee(6)(δW o;es in the notation of Rodriguez-Marek et al., 2011) andthe standard deviation of these residuals is σIII (ϕSS in thenotation of Rodriguez-Marek et al., 2011).
The adjustments in Figure 10 have been applied to the
PGAs on Figure 4 and the PGVs on Figure 5. In Figures 11
and 12, the adjusted points are compared with four NGA
models for V S30  760 m=s.
Figure 13 shows the residuals (δIIS
es ) and the event term(E2e ) for one of the NGA relations, Campbell and Bozorgnia(2008), as a function of distance. Table 4 gives, for the adjusted peak values, statistical parameters equivalent to thosefor the unadjusted residuals in Table 3.
One surprising result is that the adjustment has significantly reduced the standard deviation of the residuals.
The decrease is made possible because the adjustments
are customized to the stations. The change can be seen by

−93−

1944

J. G. Anderson, H. Kawase, G. P. Biasi, J. N. Brune, and S. Aoi
2

GM PGV (Adjusted to 760 m/s), cm/s

10

1

10

0

10

−1

10

−2

10

Adjusted Data
Abrahamson and Silva (2008)
Boore and Atkinson (2008)
Campbell and Bozorgnia (2008)
Chiou and Youngs (2008)
1

2

10

3

PGV, adjusted to V S30  760 m=s, using the adjustment factors from Kawase and Matsuo (2004a,b). The NGA model
predictions in this figure also use V S30  760 m=s.

Campbell & Bozorgnia
2
1

Data
exp(μ) = 1.18

0
−1
0

2

10

20

40
60
rflt, km

80

100

1.5

1

10

1
δ − pgv

GM PGA (Adjusted to 760 m/s), cm/s2

10

10Figure 12.

δ − pga

Figure 10. Custom station adjustments for PGA and PGV from
Kawase and Matsuo (2004a,b), used in this paper to adjust observed
values to a site with V S30  760 m=s. The adjustments are plotted asa function of estimated values of V S30 for each site, as explained in
the text. Estimates of V S30 are not used in this paper, but are used by
Atkinson and Casey (2003). Adjustment factors attributed to Atkinson and Casey (2003) are developed for response spectra, not PGA
or PGV. For PGA, (a) their adjustments for 5 Hz, and for PGV,(b) their adjustments for 2.5 Hz.

3

10

rflt, km

0

10

Data
exp(μ) = 1.46

0.5
0

−1

10

−2

10

Adjusted Data
Abrahamson and Silva (2008)
Boore and Atkinson (2008)
Campbell and Bozorgnia (2008)
Chiou and Youngs (2008)
1

−0.5
0

2

10

10

3

10

r , km
flt

20

40
60
rflt, km80

100

Figure 13. Residuals of PGA and PGV to the Campbell and Bozorgnia (2008) ground-motion prediction equation. As shown in
Table 4, this model has the lowest residuals of the four NGA relations.

Figure 11.

PGA, adjusted to V S30  760 m=s, using the adjustment factors from Kawase and Matsuo (2004a,b). The NGA model
predictions in this figure also use V S30  760 m=s.

comparing σIIE in Table 3 with σIII in Table 4. The standard
deviations decrease from ∼0:59 to ∼0:42 for PGA, and from
∼0:53 to ∼0:35 for PGV. In their pilot study of peak motions
recorded in Mexico, Anderson and Uchiyama (2011) found acomparable σIII  0:42 for PGA, and larger σIII  0:46 for
PGV. Rodriguez-Marek et al. (2011) found a larger value forPGA (0.497). Overestimates of seismic hazard can result
from the ergodic assumption (Anderson and Brune, 1999a),and examples such as this of single-station sigma provide

annecdotal examples on the lower limits of sigma for future
GMPE development.

This adjustment also has a significant effect on theevent terms. For PGA and PGV, the mean event terms havebeen decreased by a factor of ∼2:1 and increased by ∼10%
respectively. In spite of these adjustments, the event terms
are still positive. Table 4 finds that the PGAs average about
1.3 times greater than the model (∼  1:7σEvents;a ), and that
the PGVs average about 1.7 times greater than the model(∼  2:2σEvents;v ). These event terms are still high, but theestimates of median acceleration are greatly improved. It
is worth noting that the adjustment by site class of Atkinson−94−Ground Motions in the Fukushima Hamadori, Japan, Normal-Faulting Earthquake

1945

Table 4
Event Terms for the NGA GMPEs, after Adjusting for the Kawase and Matsuo (2004a)
Site Response Terms
Acceleration

Velocity

GMPE

Ee

εa

σIII

eEe

Ee

εv

σIII

e Ee

Abrahamson and Silva (2008)
σIa  0:586 σIv  0:583
σEvents;a  0:326 σEvents;v  0:319
Boore and Atkinson (2008)
σIa  0:564 σIv  0:560
σEvents;a  0:260 σEvents;v  0:256
Campbell and Bozorgnia (2008)
σIa  0:524 σIv  0:525
σEvents;a  0:219 σEvents;v  0:203
Chiou and Youngs (2008)
σIa  0:557 σIv  0:545
σEvents;a  0:271 σEvents;v  0:241
Average

0.35

0.60

0.42

1.42

0.52

0.90

0.35

1.69

0.11

0.20

0.41

1.12

0.59

1.05

0.34

1.80

0.16

0.31

0.42

1.18

0.38

0.72

0.35

1.46

0.58

1.04

0.47

1.79

0.69

1.26

0.37

2.00

1.35

0.54

0.30

and Casey (2003), which on average exceed our adjustments(Fig. 10), reduced event terms to near zero for the event intheir study.
All four of the NGA relations have a term for normalfaulting mechanisms. The terms for normal faulting modify
the predicted log (base e) of the PGA in the Abrahamson andSilva (2008), Boore and Atkinson (2008), Campbell and Bozorgnia (2008), and Chiou and Youngs (2008) models by−0:06, −0:25, −0:12, and −0:25, respectively, in which
the negative sign indicates reduction. The respective adjustments for PGV are −0:06, −0:45, 0.0, and −0:06. For PGA,
without the normal-faulting adjustments the Boore and Atkinson (2008) model would have essentially no residual, andCampbell and Bozorgnia (2008) event term would be smaller
than the interevent standard deviation, σEvents;a . For PGV,
without the normal-faulting term the Boore and Atkinson(2008) model would overestimate the adjusted data. With
these exceptions, the normal-faulting adjustments are not, bythemselves, possible explanations for underestimating the
ground motions in the Fukushima Hamadori earthquake.
Stress Drop
There are several ways that stress drop is measured.
Studies of the Fukushima Hamadori earthquake allow estimates of the static stress drop, and provide estimates andcomparisons with other earthquakes of stress drop determined from energy and high-frequency spectral amplitudes.
To evaluate the static stress drop we consider how theseismic moment might be divided between the two overlapping ruptures. With overlapping ruptures, it is unclear how toestimate the static stress drop for this event, so we consider
some alternative models. These are presented in Table 5. The
two most commonly used rupture models are the circularfault model and the dip-slip fault model (Kanamori and Anderson, 1975). The model for a dip-slip earthquake assumes
that the fault is long compared to its width. Because that is

1.72

not the case for either of the faults that ruptured in the
Fukushima Hamadori earthquake, we use the equation fora circular fault,

7π D
μ E;
Δτ S 
(7)16 RE
 E is the mean slip, andwhere Δτ S is the static stress drop, D
RE is the radius of a circular fault with the same area as therectangular approximations given in Table 2. With theseassumptions, Table 5 gives the mean slip and the stress drop
for five alternative distributions of the seismic moment from
Table 1 between the two active faults.
Because they overlap, extreme models would be to put
all the moment on one, and to consider the other to be sympathetic (Fraction  1:0 cases in Table 5). These lead toΔτ S  6:7 MPa if all the moment is released on the Idosawa
fault, and Δτ S  8:4 MPa if entirely on the Yunotake fault.
A model that divides the moment in proportion to the maximum surface slip, which was 2.0 m on the Idosawa fault and
0.8 m on the Yunotake fault, would have fslip  2:0=2:0 
0:8  0:71 of the slip would be on the Idosawa fault. Thisis close to dividing the moment 75:25 between the Idosawa
and Yunotake faults, and leads to Δτ S  5:0 MPa andΔτ S  2:1 MPa, respectively. Dividing the moment equally
between the two faults leads to Δτ S  3:3 and 4.2 MPa,
respectively. In summary, assuming both faults contributed−95−

Table 5
Mean Slip and Stress-Drop Estimates for the Two Faults Activein the Fukushima Hamadori Earthquake
Fraction of Total Moment on the Idosawa Fault
Fault

Idosawa
Yunotake E , cm
D
Δτ S , MPa
 E , cm
D
Δτ S , MPa

1.0

0.75

0.5

0.25

0.0

149
6.7

112
5.0
43
2.1

75
3.3
87
4.2

37
1.7
131
6.3

174
8.4

1946

J. G. Anderson, H. Kawase, G. P. Biasi, J. N. Brune, and S. Aoi

to the moment, the static stress drops on each fault are in therange from 2–5 MPa.
A global study of earthquakes based on teleseismic observations (Allmann and Shearer, 2009) found that the stressdrop for normal-faulting mechanisms is distributed primarily
from 1–10 MPa, with a mean value apparently about 3 MPa.
This suggests that the static stress drop on the two primary
faults in the Fukushima Hamadori earthquake is similar toother normal faults.
Another potentially useful way to compare with otherresults would be based on results of Hecker et al. (2010).
This study determines the ratio of maximum displacement torupture length. From Table 2, this ratio is 1:5 × 10−4 for the
Idosawa fault, and 0:5 × 10−4 for the Yunotake fault. The
former is in a range of ratios typical of faults with smaller
cumulative displacement in Hecker et al. (2010), whereas
the latter is lower than most faults in their data set. In summary, the static stress drop associated with the Fukushima
Hamadori earthquake is in a range that is typical of otherearthquakes with normal-faulting mechanisms.
Nakano (2013) determined stress drop for several of
the earthquakes in the Fukushima Hamadori sequence. Extending the approach of Kawase and Matsuo (2004a), they
combined the seismic moment determined from regional
recordings based on F-NET data (Okada et al., 2004) and corner frequency measured from local stations to find the stressdrop using the equations in Brune (1970). By their results,
ΔσBrune  3:1 MPa. This is within the range of the stressparameters found by Nakano (2013) for other Japanese earthquakes with Mw ≳ 5:3, with reverse or strike-slip mechanisms.
Satoh and Tsutsumi (2012) measured the short-period
spectral level of the Hamadori earthquakes and comparedthat with other Japanese events. According to their model,
the short-period level is proportional to the product of astress-drop parameter (not necessarily the same as ΔσBrune )and fault radius. They find that their stress parameter is approximately constant for events with Mw ≳ 5:7, with reversefaulting earthquakes having a higher stress parameter, on
average, than strike-slip-faulting earthquakes. They find that
the Fukushima Hamadori mainshock has a stress parameter,
which is intermediate between the overlapping populationsof reverse and strike-slip stress parameters.
Finally, Somei et al. (2012) find stress drop using thespectra of coda waves to find seismic moment and cornerfrequency, following the approach of Mayeda and Walter(1996). They find relatively small scatter in stress-drop estimates for M w ≳ 4:5. In this range, and by this method, theyfind that the stress drop of strike-slip- and reverse-mechanism
earthquakes overlap almost completely, and that the Fukushima Hamadori normal faulting event is near the center ofthat range.
In summary, it appears that the static stress drop of the
Fukushima Hamadori earthquake is consistent with theglobal population of normal-faulting earthquakes and thatthe other stress parameters are consistent with the Japanese
population of strike-slip and reverse earthquakes.

Simultaneous Rupture of the Two Main Faults
The results presented in the previous section indicatethat the static stress drop of the two faults, taken individually,is fairly typical of normal-faulting earthquakes. However,given the overlap of the geometry, if the two faults ruptured
simultaneously the combined radiation from the two ruptures
might be a cause of the above-average ground motions.
Hikima (2012) inverted the strong-motion waveforms tofind a kinematic model of slip in the earthquake. He findsthat the Idosawa fault slipped first, and the Yunotake fault
slipped second. In this model, the Idosawa fault slip takes
place in the first 8 seconds after the origin time, with thelargest slip in the time interval from 3–7 seconds. The major
rupture on the Yunotake segment takes place in the intervalfrom about 9–13 seconds.
Shiba and Noguchi (2012) also found a kinematic modelfor the source based on waveform inversion. In their model,
slip also occurs first on the Idosawa fault, and has largely
ended on that fault after about 6 s. Rupture on the Yunotake
fault in this model starts at about 5 s and ends at about 13 s
after the origin. Thus, their model allows simultaneity of rupture during only a small part of the earthquake.
Both of these studies, then, find sequential rupture of
the two segments, rather than simultaneous rupture. Given
the small separation of the two faults, there would be little
overlap in ground motion from the two segments, probably
not significantly more than the overlap in ground motion
from multisegment ruptures of faults in which subsequent
segments follow a more linear geometry.
GMPE Calibrations
The fourth hypothesis for the high event terms that needs
to be considered closely is the calibration of the GMPEs. As
noted in Figure 1, the data at large magnitude and small distance used to calibrate the GMPEs for the effect of normal
faulting in the NGA relations is sparse. In the preliminary
findings above, other obvious factors do not seem to explain
the somewhat high event terms, leaving the calibration of the
GMPEs as an obvious potential contributor.
Additional Considerations
There are two themes of evidence that might support theidea that the NGA relations are reasonably correct and theevent terms are high for this particular event.
Precarious Rocks. Brune (2000, 2003) discusses precariously balanced rocks (PBRs) in normal-faulting and transtensional environments. For the sake of this discussion, a
PBR is a rock that is balanced on, but not mechanically attached to, its pedestal, and has a high enough ratio of height
to width that it is easily toppled by strong ground motions(Anderson et al., 2011). For instance, some of the PBRs thathave been cited to constrain ground motions 15 km from the
San Andreas fault are expected to topple in ground motions−96−

Ground Motions in the Fukushima Hamadori, Japan, Normal-Faulting Earthquakewith PGAs exceeding ∼0:4g. Brune (2000) observed equallyfragile PBRs located on the footwall close to several normal
faults with Quaternary activity. Based on this, and supportedby the physical model results of Brune and Anooshehpoor(1999), he proposed the hypothesis that normal-faulting
events at short distances on the footwall (e.g., less than 3 km
from the fault trace) do not, in general, cause ground motions
sufficient to topple many precarious rocks, as they would at acomparable distance from the trace of a strike-slip or thrustfault. Brune (2003) suggests that ground motions also tend tobe low near transtensional stepovers of strike-slip faults. Of
the three examples presented there, the Honey Lake zone ofprecarious rocks is most relevant (Briggs et al. [2013]). The
rocks there are on the footwall of a normal fault associatedwith a major extensional stepover of strike-slip faults.
Unfortunately, none of the ground motions observedin the Fukushima Hamadori earthquake are located in a position to test the hypothesis in Brune (2000), and this earthquake has the wrong mechanism to test the Brune (2003)
hypothesis. Only two K-NET stations (FKS010 and
FKS011; Fig. 8) and one KiK-net station (FKSH14; Fig. 6)
are at an azimuth, which can be considered to sample thefootwall ground motions, and all three are over 10 km fromthe vertical projection of the fault. The K-NET stations bothhave relatively high values of V S30 , and the KiK-net station
has a relatively low surface/downhole ratio. Thus, all three
are expected to have somewhat lower ground motions basedon site condition properties. To the extent that these threestations have below-average ground motions, the cause isambiguous at best and would require more study.
A loose extrapolation of the Brune (2000, 2003) studies,
however, leads one to consider the hypothesis that groundmotions in normal and transtensional environments are generally lower than in strike-slip environments. This hypothesis
is consistent with the normal-faulting adjustments that werediscussed previously. Thus, it is reasonable to ask whether
the observed ground motions imply that precarious rocks
any place would be toppled. Precarious rocks are situated onsolid rock pedestals, so correction of ground motions torock site conditions (V S30  760 m=s) is essential to compare with their hypothesis.
All but one of the PGAs are corrected to below 0:3g
(Fig. 11). The lone exception is station FKSH13. This station
is above the deepest part of the fault on the hanging wall(Fig. 6), and thus does not test the Brune (2000) hypothesis.
We consider whether ground motion on rock at this location
would have caused precarious rocks to be toppled. The initialadjustment for site response would seem to suggest groundmotions capable of causing at least some precarious rocks totopple. The geometric mean horizontal PGA at site FKSH13
is adjusted from 922 cm=s2 to 709 cm=s2 . The geometricmean horizontal PGV at site FKSH13 is practically unchanged, adjusting from 30:6 cm=s to 30:5 cm=s2 . The toppling of precarious rocks depends on both parameters. The
value of PGA is essential as a threshold to initiate rocking anddepends on the shape of the rock. The additional PGV thresh-

1947

Figure 14.

Surface and downhole accelerogram from surface
sensors at FKSH13. The vertical component of the downhole sensor
did not produce a seismogram.

Figure 15. Surface and downhole velocity at FKSH13. The
vertical downhole sensor did not produce a seismogram.
old (or a related measure of longer period motions such as
spectral acceleration for an oscillator period of 1.0 s) depends
on the size of the rock. A PGV of 30 cm=s is sufficient tocause smaller precarious rocks to topple, but large rocks(e.g., more than 2 m high) may require higher velocities.
Figure 14 shows the surface and downhole acceleration,
Figure 15 shows the surface and downhole velocity, andFigure 16a shows the Fourier spectra of the record from station FKSH13. The vector horizontal PGA at the surface was1306 cm=s2 , whereas downhole it was 179 cm=s2 (Fig. 16b).
The Fourier spectra show a very broad peak of amplificationfrom 3–12 Hz, peaking at 6–7 Hz, that is responsible for thelarge surface acceleration. The velocity seismogram shows
prominent high-frequency energy at this spectral peak, similar to IBRH12 in Figure 7. Thus, this location appears tohave a contribution to shaking from resonant site conditions.
A station on solid rock might have expected acceleration
two–three times the downhole value, but the peak may be

−97−

1948

J. G. Anderson, H. Kawase, G. P. Biasi, J. N. Brune, and S. Aoi

Figure 17.
Figure 16.

(a) Smoothed horizontal surface and downhole
spectra at FKSH13. (b) Ratio of horizontal surface to horizontal
downhole Fourier spectra.

Cabin adjacent to fault that ruptured in the 1954
Dixie Valley earthquake (M w 6.9). The photograph was taken on
20 December 1954 by Karl Steinbrugge. The man in the photois Vincent Bush, a member of the team investigating the earthquake.
Source: Karl V. Steinbrugge Collection: S653, National Information
Service for Earthquake Engineering (NISEE), University of California, Berkeley. The color version of this figure is available only in theelectronic edition.

smaller than the adjusted PGA based on the Kawase andMatsuo (2004a,b) correction. A correction of the seismogram in the frequency domain would be needed to more accurately predict these effects, but such a correction is beyond
the scope of this paper. All of these things considered, asolid rock outcrop at the location of this station would likely
have accelerations capable of toppling most precarious
rocks. Thus, an extrapolation of the Brune (2000, 2003) hypotheses to infer reduced shaking at all distances from a normal fault would not be supported by this data.
1954 Cabin. Finally, we consider Figure 17, a classic picture of a normal-fault scarp in Nevada, formed in the M w 6.9
Dixie Valley earthquake of 16 December 1954. Figure 18
shows a different view of the cabin, which shows that thecabin is in a graben caused by the faulting. A remarkable
observation by the scientists investigating the faulting wasthat except for the tilt due to differential settlement near thefault, the structure was essentially undamaged. Furthermore,inside the structure a tin cup had not even been shaken off the
shelf (Slemmons, 1957; Steinbrugge and Moran, 1957). As
with the precarious rocks discussed by Brune (2000), thisevidence suggests that motions near the fault were not particularly severe, even as the ground failed on the hanging
wall. Slemmons (1957) also noted that in the Dixie Valley
settlement three miles (∼5 km) east of this cabin, above thehanging wall, “the up-and-down motion was so strong as tothrow some of the residents out of bed and then repeatedly
knock them to the floor as they made their way outdoors.” Healso notes that the movement of objects inside buildings was
spatially inconsistent. At the Dixie Valley settlement, on
sediments, Slemmons reports at least one home in which
essentially everything that could be moved was shaken to thefloor, as would be expected from an acceleration of thenature of what was observed at FKSH13.

Figure 18. Karl Steinbrugge standing next to the cabin at the1954 Dixie Valley rupture (M w 6.9). The photograph was taken on
20 December 1954 by Donald Moran. Source: Karl V. Steinbrugge
Collection: S662, National Information Service for Earthquake Engineering (NISEE), University of California, Berkeley. The color
version of this figure is available only in the electronic edition.
Accelerations and velocities that are much stronger
above the hanging wall are predicted by a lattice model ofnormal faulting (Shi et al., 2003). The lattice model also predicts that these strong motions extend all the way to the fault
scarp. Thus, the inferred low amplitudes very near the fault
trace at the cabin may be a local effect brought about by
highly nonlinear response of the sediments in the fault zone
(as can be observed by the shattered surface around the cabinin Fig. 17), whereas motions farther away may be larger.
Indeed, besides this anecdotal evidence of strong shaking on the hanging wall a few kilometers from the fault trace,
the lack of precarious rocks provides some support for thesuggestion that low accelerations do not persist to larger−98−Ground Motions in the Fukushima Hamadori, Japan, Normal-Faulting Earthquake
distances on the hanging wall. Purvance et al. (2008) find precarious rocks that constrain past ground motions in a granite
outcrop on the hanging wall of one splay of the active Genoa
fault in western Nevada. However, one of us (JNB) has foundthat similar granites located 5–10 km from the fault on the
hanging wall have no precarious rocks, and indeed rather havethe appearance of extensive rockfalls such as might be expected to be triggered by strong shaking. All of this evidence
is consistent with the idea that the frictional forces, which must
be overcome at depth to cause a normal-faulting earthquake,and the consequent seismic radiation, are comparable to theforces associated with other mechanisms.
We are not aware of any detailed intensity surveys near
the faults that broke in the Fukushima Hamadori earthquake.
Koarai and Okaya (2011) and Hisada (2012) reported a small
number of Japanese wooden houses, which were heavily
damaged, but almost all of them seem to be associated witheither fault movement or ground failure (including liquefaction) so that they are not indicative of severe ground motions.
At first glance, this lack of heavy damage to wooden houses
elsewhere may suggest weak ground motions, contrary to the
Dixie Valley observation. However, to create heavy damage
to Japanese wooden houses we need both high PGA andhigh PGV, namely PGA ≥ 0:8g and PGV ≥ 100 cm=s basedon the damage survey in Kobe (Kawase, 2011). From theKiK-net and K-NET records (Figs. 6 and 8) PGV in the source
region is not observed in excess of 53 cm=s, which is well
below 100 cm=s, so the absence of severe damage to wooden
houses is exactly what we should expect. So far as we know, adetailed survey of the effect of the shaking on contents of thestructures in the region has not been carried out yet, so a basis
for comparison with Dixie Valley apparently does not exist.

Summary
To sample two parts of the seismic spectrum in a preliminary assessment of the potential impact of this eventon prediction of median estimates of ground motion innormal-faulting earthquakes, this paper compared observed
PGA and PGV with the median values of four NGA relations.
The models for PGA are generally lower than the observations in the Fukushima Hamadori earthquake by a factorof ∼2:3–3:7, depending on the GMPE, at all distances to
100 km. For PGV the differences are smaller, with the data
exceeding the predictions by factors of 1.4–1.8.
We consider four hypotheses to explain this. The firstis that the Fukushima Hamadori earthquake is unusual fornormal faulting. The case for this is perhaps easily made;it is triggered by a large subduction-zone earthquake, so it
might be easy to presume it should be different from otherevents. However, aside from this unusual triggering, we areunable to identify physical parameters describing the faultingthat set this earthquake apart. Several stress-drop measurements of the Fukushima Hamadori earthquake are similar tomeasurements from other crustal earthquakes in Japan. The
comparison events predominantly have strike-slip or reverse

1949

mechanism. Second, we find that the available kinematic
models of the Fukushima Hamadori earthquake show thetwo participating faults ruptured sequentially, so simultaneous contributions from those two faults does not appear to bea viable explanation for the high event terms.
Another alternative is that the records of the event havehigher motions than average due to a prevalence of sites withstrong resonances. Our corrections are based on customized
station terms found by Kawase and Matsuo (2004a,b). Thesestation terms cause a significant reduction in the standard
deviation of corrected data from the GMPEs, and thus provide
strong support for the idea that corrections based on datafrom smaller earthquakes can improve ground-motion prediction over models based on simplified site parameters suchas V S30 . These corrections also strongly reduced the event
term for PGA, indicating the power of customized site terms
to impact mean ground-motion predictions. The event terms
remaining after this correction are sensitive to the empiricaladjustment of Kawase and Matsuo between V S30  3450 m=s
and V S30  760 m=s. Because the data used by Kawase andMatsuo (2004a,b) are scattered, an objective for future studies
could be to seek additional theoretical guidance for the formand size of that adjustment, and to understand the differencesfrom adjustments used by Atkinson and Casey (2003).
More sophisticated site correction is still possible, suchas adjustment of records for Fourier spectral amplifications
determined by joint inversions for source, path, and site effects (e.g., Kawase, 2006). Thus, the test of the effect of site
response in this paper is incomplete, but if we accept the preliminary results we would conclude that site effects accountfor a large part of the large event term for PGA. If we presume
that adjustments for site resonance could completely explain
the high event term, it raises complex issues for ground-motion prediction. It casts doubt on the efficacy of V S30 alone tomodel these effects even to first order, and implies resonanceneeds to be built into the models. To achieve that, the model
will need to be based on perterbations from a reference-site
velocity and density profile.
The last alternative is that GMPEs have underestimatedaccelerations for large normal-faulting events at distanceslarger than about 5 km, for which the instrumental data fromthe Fukushima Hamadori earthquake provide control. Observational geologic data, mainly from precariously balanced
rocks, suggest that this may not be the case at closer distances. An underestimate of ground motions at high frequencies
would primarily affect hazard estimates on the U.S. National
Seismic Hazard Map in regions dominated by normal faulting, such as the Basin and Range region. To the extent thatthese Japanese data influence future GMPEs, this might resultin a small increase in the hazard in a future version of the
U.S. national map.

Data and Resources
The NIED website for site velocity profiles and accelerograms from K-NET and KiK-net stations is now

−99−

1950

J. G. Anderson, H. Kawase, G. P. Biasi, J. N. Brune, and S. Aoi(http://www.kyoshin.bosai.go.jp/), last accessed November
2012. We accessed the same data from predecessors to the current site: for K-NET: (http://www.k-net.bosai.go.jp/, last accessed January 2012) and for KiK-net (http://www.kik.bosai.go.jp/kik/index_en.shtml, last accessed January 2012).
GMPEs of Abrahamson and Silva (2008), Boore and Atkinson(2008), Chiou and Youngs (2008), and Campbell and Bozorgnia (2008) were calculated using the Attenuation Relationship
Plotter in OPEN_SHA (Field et al., 2003; http://www.opensha.org/, last accessed January 2012). The NGA flatfile is available
from PEER (http://peer.berkeley.edu/nga/flatfile.html, last accessed August 2012). Strong-motion data from Italy is available from http://itaca.mi.ingv.it/ItacaNet/itaca10_links.htm(last accessed January 2013), and from Greece is available
from http://www.itsak.gr/en/page/data/strong_motion/ (last accessed January 2013). The Global Centroid Moment Tensor
Project database was searched using www.globalcmt.org/
CMTsearch.html (last accessed January 2013).

Boore, D. M., and G. M. Atkinson (2008). Ground-motion prediction
equations for the average horizontal component of PGA, PGV, and
5%-damped PSA at spectral periods between 0.01 s and 10 s, Earthq.
Spectra 24, 99–138.
Boore, D. M., and W. B. Joyner (1997). Site amplifications for generic rocksites, Bull. Seism. Soc. Am. 87, 327–341.
Boore, D. M., E. M. Thompson, and H. Cadet (2011). Regional correlationsof V S30 and velocities averaged over depths less than and greater than
30 meters, Bull. Seismol. Soc. Am. 101, 3046–3059.
Briggs, R. W., S. G. Wesnousky, J. N. Brune, M. D. Purvance, and S. A.
Mahan (2013). Low footwall accelerations and variable surface rupture
behavior on the Fort Sage Mountains fault, northeast California, Bull.
Seism. Soc. Am. 103, no. 1, 157–168.
Brune, J. N. (1970). Tectonic stress and the spectra of seismic shear waves
from earthquake, J. Geophys. Res. 75, 4997–5009.
Brune, J. N. (2000). Precarious rock evidence for low ground shaking on
the footwall of major normal faults, Bull. Seismol. Soc. Am. 90,
1107–1112.
Brune, J. N. (2003). Precarious rock evidence for low near-source accelerations for transtensional strike-slip earthquakes, Phys. Earth Planet.
In. 137, 229–239.
Brune, J. N., and A. Anooshehpoor (1999). Dynamic geometrical effects onstrong ground motion in a normal fault model, J. Geophys. Res. 104,
809–815.
Campbell, K. W., and Y. Bozorgnia (2008). NGA ground motion model forthe geometric mean horizontal component of PGA, PGV, PGD and 5%
damped linear elastic response spectra for periods ranging from 0.01 to
10 s, Earthq. Spectra 24, 137–171.
Chiou, B. S.-J., and R. R. Youngs (2008). An NGA model for the average
horizontal component of peak ground motion and response spectra,
Earthq. Spectra 24, 173–215.
Chiou, B., R. Darragh, N. Gregor, and W. Silva (2008). NGA Projectstrong-motion database, Earthq. Spectra 24, 23–44.
Field, E. H., T. H. Jordan, and C. A. Cornell (2003). OpenSHA: A developing community-modeling environment for seismic hazard analysis,
Seismol. Res. Lett. 74, no. 4, 406–419.
Geology Research Team (2011). Landforms, geology, and active faults (inJapanese), Powerpoint presentation, download from http://www.pwri.go.jp/eng/news/2011/0602/ (last accessed 4 May 2012).
Hecker, S., T. E. Dawson, and D. P. Schwartz (2010). Normal-faulting
slip maxima and stress-drop variability, Bull. Seismol. Soc. Am. 100,
3130–3147.
Hikima, K. (2012). Rupture process of the April 11, 2011 Fukushima
Hamadori earthquake (M j 7.0), Zishin 64, 243–256 (in Japanese).
Hisada, Y., T. Kubo, K. Matsuzawa, T. Matsumoto, A. Tanabe, and A.
Morikawa (2012). Investigation of building damage near the surface
faults of 2011 Fukushima-ken-Hama-Dori earthquake, special issue, J.
Japan Assoc. Earthq. Eng. 12, no. 4, 104–126, doi: 10.5610/
jaee.12.4_104.
Imanishi, K., R. Ando, and Y. Kuwahara (2012). Unusual shallow normalfaulting earthquake sequence in compressional northeast Japan activated after the 2011 off the Pacific coast of Tohoku earthquake,
Geophys. Res. Lett. 39, L09306, doi: 10.1029/2012GL051491.
Kanamori, H., and D. L. Anderson (1975). Theoretical basis of some empirical relations in seismology, Bull. Seismol. Soc. Am. 65, 1073–1095.
Kato, A., S. Sakai, and K. Obara (2011). A normal-faulting seismic sequence
triggered by the 2011 off the Pacific coast of Tohoku earthquake:
Wholesale stress regime changes in the upper plate, Earth Planets
Space 63, no. 7, 745–748, doi: 10.5047/eps.2011.06.014.
Kawase, H. (2006). Site effects derived from spectral inversion method for
K-NET, KiK-net, and JMA strong-motion network with special reference to soil nonlinearity in high PGA records, Bull. Earthq. Res. Inst.
Tokyo Univ. 81, 309–315.
Kawase, H. (2011). Strong motion characteristics and their damage impact
to structures during the Off Pacific Coast of Tohoku Earthquake ofMarch 11, 2011: How extraordinary was this M 9.0 earthquake?, inProc. of the 4th IASPEI/IAEE International Symposium on Effects

Acknowledgments
The authors thank Hiroe Miyake and Hongjun Si for providing information, which was not otherwise available, and for several very helpful interactions. John Anderson thanks the Earthquake Research Institute of the
University of Tokyo for supporting a sabbatical visit and a return visit toTokyo, which provided invaluable insights and interactions making thispaper possible. Burt Slemmons identified critical references for the 1954
earthquake, and helped identify the man in the cabin-fault photo. Adrian
Rodriguez-Marek provided helpful comments on an early version of thispaper. This paper also benefited from helpful reviews by Gail Atkinson, Ivan
Wong, and an anonymous reviewer. We thank the National Research Institute for Earth Science and Disaster Prevention for the making the K-NET and
KiK-net strong-motion data and downhole profiles freely available. This
study was completed with partial support from the U.S. Geologial Survey
NEHRP External Grants Program, under Award Number G12AP20024.

References
Abrahamson, N. A., and W. Silva (2008). Summary of the Abrahamson &Silva NGA ground-motion relations, Earthq. Spectra 24, 67–97.
Aki, K., and P. G. Richards (2002). Quantitative Seismology, Second Ed.,
University Science Books, Sausalito, California, 700 pp.
Allmann, B. P., and P. M. Shearer (2009). Global variations of stress drop formoderate to large earthquakes, J. Geophys. Res. 114, B01310, doi:
10.1029/2008JB005821.
Anderson, J. G. (2013). Surface motions on near-distance rock sites in the2011 Tohoku earthquake, Earthq. Spectra 29, no. S1, S23–S35.
Anderson, J. G., and J. N. Brune (1999a). Probabilistic seismic hazard analysis without the ergodic assumption, Seismol. Res. Lett. 70, 19–28.
Anderson, J. G., and Y. Uchiyama (2011). A methodology to improve groundmotion prediction equations by including path corrections, Bull. Seismol. Soc. Am. 101, no. 4, 1822–1846, doi: 10.1785/0120090359.
Anderson, J. G., J. N. Brune, G. Biasi, A. Anooshehpoor, and M. Purvance(2011). Workshop report: Applications of precarious rocks and related
fragile geological features to U. S. National Hazard Maps, Seismol.
Res. Lett. 82, no. 3, 431–441, doi: 10.1785/gssrl.82.3.431.
Atkinson, G. M., and R. Casey (2003). A comparison of ground motions
from the 2001 M 6.8 in-slab earthquakes in Cascadia and Japan, Bull.
Seismol. Soc. Am. 93, 1823–1831.
Bindi, D., F. Pacor, L. Luzi, R. Puglia, M. Massa, G. Ameri, and R. Paolucci(2011). Ground motion prediction equations derived from the Italianstrong motion database, Bull. Earthq. Eng. 9, 1899–1920, doi:
10.1007/s10518-011-9313-z.

−100−Ground Motions in the Fukushima Hamadori, Japan, Normal-Faulting Earthquakeof Surface Geology on Seismic Motion, August 23–26, 2011,
University of California Santa Barbara, Santa Barbara, California.
Kawase, H., and H. Matsuo (2004a). Amplification characteristics of
K-NET, KiK-NET, and JMA Shindokei network sites based on the
spectral inversion technique, in 13th World Conference on Earthquake
Engineering, Vancouver, Canada, Paper No. 454.
Kawase, H., and H. Matsuo (2004b). Relationship of S-wave velocity structures and site effects separated from the observed strong motion dataof K-NET, KiK-net, and JMA network, J. Japan Assoc. Earthq. Eng. 4,
no. 4, 126–145 (in Japanese).
Kinoshita, S. (1998). Kyoshin net (K_NET), Seismol. Res. Lett. 69, 309–332.
Kinoshita, S. (2005). Development of strong-motion observation network
constructed by NIED, in Directions in Strong Motion Instrumentation,
P. Gulkan and J. G. Anderson (Editors), NATO Science Series IV,
Earth and Environmental Sciences, Vol. 58, Springer, Dordrecht,
the Netherlands, 181–196.
Koarai, M., and T. Okaya (2011). Field Survey Report (Supplement) on the
Seismic Damage around Iwaki City (in Japanese), downloaded fromhttp://www.gsi.go.jp/common/000060960.pdf (last accessed 24 June2012).
Luzi, L., F. Sabetta, F. Mele, and B. Castello (2010). Italian strong motion
database relative to the period 1972–2004: Motivations and aims, Bull.
Earthq. Eng. 8, 1159–1174, doi: 10.1007/s10518-009-9140-7.
Maruyama, T., A. Saito, T. Azuma, K. Taniguchi, M. Yoshimi, and T.
Hayashida (2011). Reconnaissance Survey Report on the Surface Ruptures along the Idosawa Fault Associated with the Fukushima Hamadori Earthquake of April, 11, 2011 (in Japanese), downloaded fromhttp://unit.aist.go.jp/actfault‑eq/Tohoku/report/idosawa/idosawa.pdf(last accessed 4 May 2012).
Massa, M., F. Pacor, L. Luzi, D. Bindi, G. Milana, F. Sabetta, A. Gorini, and
S. Marcucci (2010). The ITalian ACcelerometric Archive (ITACA):
Processing of strong-motion data, Bull. Earthq. Eng. 8, 1175–1187,
doi: 10.1007/s10518-009-9152-3.
Mayeda, K., and W. R. Walter (1996). Moment, energy, stress drop, and source
spectra of western U.S. earthquakes, J. Geophys. Res. 101, 11,195–11,208.
Moya, A., and K. Irikura (2003). Determination of site effects and Q factor
using a reference event, Bull. Seismol. Soc. Am. 93, 1730–1745.
Nakano, K. (2013). Statistical properties of strong ground motions based on
the spectral inversion method, Master’s Thesis, 28 February 2013,
Dept. of Engineering, Kyoto University, Kyoto, Japan (in Japanese).
Okada, Y., K. Kasahara, S. Hori, K. Obara, S. Sekiguchi, H. Fujiwara, andA. Yamamoto (2004). Recent progress of seismic observation networks in Japan, Hi-net, F-net, K-NET and KiK-net, Earth Planets
Space 56, xv–xxviii.
Pacor, F., R. Paolucci, L. Luzi, F. Sabetta, A. Spinelli, A. Gorini,
M. Nicoletti, S. Marcucci, L. Filippi, and M. Dolce (2011). Overview
of the Italian strong motion database ITACA 1.0, Bull. Earthq. Eng. 9,
1723–1739, doi: 10.1007/s10518-011-9327-6.
Petersen, M. D., A. D. Frankel, S. C. Harmsen, C. S. Mueller, K. M. Haller, R.
L. Wheeler, R. L. Wesson, Y. Zeng, O. S. Boyd, D. M. Perkins, N. Luco,E. H. Field, C. J. Wills, and K. S. Rukstales (2008). Documentation forthe 2008 Update of the United States National Seismic Hazard Maps,
U.S. Geological Survey Open-File Report 2008–1128, 61 pp.
Power, M., B. Chiou, N. Abrahamson, Y. Bozorgnia, T. Shantz, and C. Roblee(2008). An overview of the NGA project, Earthq. Spectra 24, 3–21.
Purvance, M. D., J. N. Brune, N. A. Abrahamson, and J. G. Anderson (2008).
Consistency of precariously balanced rocks with probabilistic seismic
hazard estimates in southern California, Bull. Seismol. Soc. Am. 98,
2629–2640.
Rodriguez-Marek, A., G. A. Montalva, F. Cotton, and F. Bonilla (2011).
Analysis of single-station standard deviation using the KiK-net data, Bull.
Seismol. Soc. Am. 101, no. 3, 1242–1258, doi: 10.1785/0120100252.
Ryder, I., A. Rietbrock, K. Kelson, R. Burgmann, M. Floyd, A. Socquet,
C. Vigny, and D. Carrizo (2012). Large extensional aftershocks inthe continental forearc triggered by the 2010 Maule earthquake, Chile,
Geophys. J. Int. 188, 879–890.

1951

Satoh, T., and H. Tsutsumi (2012). Short-period spectral level, path and siteamplification derived from strong motion records during the normalfaulting earthquakes in Hamadori, Fukushima Prefecture in 2011, J.
Japan Assoc. Earthq. Eng. 12, no. 7, 1–18 (in Japanese).
Shi, B., J. N. Brune, Y. Zeng, and A. Anooshehpoor (2003). Dynamics ofearthquake normal faulting: Two-dimensional lattice particle model,
Bull. Seismol. Soc. Am. 93, 1179–1197.
Shiba, Y., and S. Noguchi (2012). Statistical characteristics of seismic source
parameters controlling broadband strong ground motions—Investigation
based on source inversion analysis, CRIEPI Rept. N11054, 28 pp. (inJapanese).
Si, H., and S. Midorikawa (1999). New attenuation relationships for peakground acceleration and velocity considering effects of fault typeand site condition, J. Struct. Constr. Eng., Architectural Institute ofJapan, no. 523, 63–70 (in Japanese).
Si, H., and S. Midorikawa (2000). New attenuation relations for peak groundacceleration and velocity considering effects of fault type and site condition, in Proceedings, 12th World Conference on Earthquake Engineering, 30 January–4 February 2000, Auckland, New Zealand, New
Zealand Society for Earthquake Engineering, Upper Hutt, New
Zealand, Paper 0532.
Si, H., T. Kawasato, M. Ohba, T. Ishikawa, T. Masatsuki, and A. Sawaii (2012).
A comparison of existing ground motion prediction equation to normal
faulting earthquake, in Summaries of Technical Papers of Annual Meeting, Architectural Institute of Japan, B-2, 177–178 (in Japanese).
Slemmons, D. B. (1957). Geological effects of the Dixie Valley-Fairview
Peak, Nevada, earthquakes of December 16, 1954, Bull. Seism. Soc.
Am. 47, 353–375.
Somei, K., K. Asano, T. Iwata, and K. Miyakoshi (2012). Comparisons ofsource characteristics between recent inland crustal earthquake sequences
inside and outside of Niigata-Kobe tectonic zone, Japan, Poster S53A2488, American Geophysical Union Fall Annual Meeting, San Francisco.
Spudich, P., W. B. Joyner, A. G. Lindh, D. M. Boore, B. M. Margaris, andJ. B. Fletcher (1999). SEA99: A revised ground motion prediction relation
for use in extensional regimes, Bull. Seismol. Soc. Am. 89, 1156–1170.
Steinbrugge, K. V., and D. F. Moran (1957). Engineering aspects of the Dixie
Valley-Fairview Peak earthquakes, Bull. Seism. Soc. Am. 47, 335–348.
Theodulidis, N., I. Kalogeras, C. Papazachos, V. Karastathis, B. Margaris,
Ch. Papaioannou, and A. Skarlatoudis (2004). HEAD 1.0: A unifiedHEllenic accelerogram database, Seismol. Res. Lett. 75, 36–45.
Wu, C., Z. Peng, and D. Assimaki (2009). Temporal changes in site response
associated with the strong ground motion of the 2004 Mw 6.6 MidNiigata earthquake sequences in Japan, Bull. Seismol. Soc. Am. 99,
no. 6, 3487–3495, doi: 10.1785/0120090108.
Nevada Seismological Laboratory
University of Nevada
Reno, Nevada 89557
jga@unr.edu(J.G.A., G.P.B., J.N.B.)
DPRI
Kyoto University
Gokasho, Uji
Kyoto 611-0011, Japan
kawase@zeisei.dpri.kyoto‑u.ac.jp(H.K.)

Natl Res Inst Earth Sci & Disaster Prevention (NIED)
3-1 Tennodai, Tsukuba
Ibaraki 305-0006, Japan
aoi@bosai.go.jp(S.A.)−101−

Manuscript received 8 August 2012

The Composite Source Model for Broadband
Simulations of Strong Ground Motions
by John G. Anderson
INTRODUCTION

the offset on the fault as a function of location and time, andcalculations of the Green’ s function.

My goal is to generate synthetic strong ground motions that are
sufficiently realistic to be useful for engineering applications.
The composite source model (CSM) uses a kinematic source
model for rupture on a finite fault. This source is propagatedto the station using a flat-layered velocity model, scattering,
and attenuation that can be measured from independent seismological observations. A key objective is to reproduce thewave propagation within the constraints of the measured velocity and Q structure. This article summarizes the CSM as implemented for the Southern California Earthquake Center(SCEC) Broadband Platform validation exercise (Goulet et al.,2015).

HISTORY
The CSM was first described by Yu (1994) and Zeng et al.
(1994). Motivated by Frankel (1991), the model describes
the source slip function as a superposition of overlapping circular subevents of random sizes, located randomly on the fault.
Multiple source realizations naturally result in different seismograms. Zeng and Anderson (1996) and Zeng and Chen (2001)
demonstrated that it is possible to find a specific realization of
the CSM that reproduces the low-frequency accelerograms in
Northridge and Chi-Chi earthquakes and that the slip functions resemble the slip function determined by other methods.
Additional applications of the model include Khattri et al.
(1994), Yu et al. (1995), Anderson and Yu (1996), and Su,
Anderson, Ni, et al. (1998), Su, Anderson, Zeng (1998), andHartzell et al. (2011).
Variations of the original model have been proposed (e.g.,Zeng et al., 1991, 1995). For simplicity, the current implementation is close to the basic version, although some aspects of the
Broadband platform implementation have not been described
elsewhere. For this reason, this article gives a brief summary of
the model on the Broadband Platform.

THEORETICAL BASIS
Representation Theorem
The theoretical basis for understanding ground motion froma finite fault is the representation theorem (e.g., Aki and Richards, 2002). Through the representation theorem, the problem of predicting ground motions is reduced to specification of
68

Seismological Research Letters

Volume 86, Number 1

Composite Source Model Description
The slip function on the fault is treated as a superposition ofsubevents that are located at random on the fault. Each subevent can be visualized as having a circular cross section, where
subevents can overlap. The rupture time of each subevent isdetermined by its distance from the hypocenter, using a constant rupture velocity, V r . The moment of each subevent depends on the subevent radius Ri and the subevent stress dropΔτS . V r and ΔτS are model parameters, whereas the radii Ri arechosen at random from a probability distribution. The moment of the ith subevent, with radius Ri , is, from Kanamori
and Anderson (1975)
M 0i 

16
ΔτS R3i :7

1

The time function of slip on each subevent is the source
time function given by Brune (1970, 1971). To be specific,
for a subevent with radius Ri the corner frequency of the subevent is
f ic 

2:34β
;
2πRi

2

and the time function is then given by_ 0i  2πf ci 2 M 0i H ττe−2πf c τ ;M

3

in which τ  t − t i , in which t i is the start time of the subevent: t i  xhyp =vr in which xhyp is the path distance of thesubevent center to the hypocenter measured along the fault,
and vr is the rupture velocity. Zeng et al. (1991) implemented
the Sato and Hirasawa (1973) pulse instead of the Brune pulse.
The Sato and Hirasawa pulse is an option controlled by input
on the Broadband Platform, but in this article all results use theBrune pulse.
The distribution of subevent sizes, and thus of rise times
among the subevents, motivated by Frankel (1991) is controlled
by fractal dimension D. Specifically, the distribution of subeventswith radius R or greater is given by Zeng et al. (1994) as
January/February 2015

−103−

doi: 10.1785/0220140098

P  −DR − R−D
max ;
D

4

(a) 0

in which P is determined by the constraint that the sum of the
moments of the subevents equals M 0 , the moment of the earthquake being modeled. Consequently,

5

N R 

P

7M 0
3−D
3−D D ≠ 3:
16ΔτS R3−D
max − Rmin 

10

5

My calculations use D  2, because then the distribution of subevent sizes follows a Gutenberg–Richter distribution with the bvalue b  D=2  1. The parameter Rmax in equation (4) is the
maximum subevent radius. My practice is to set Rmax equal tohalf of the smaller fault dimension. The minimum subevent radius Rmin is currently chosen by Rmin  Rmax =20. The momentof a subevent with radius Rmax is 8000 times greater than a subevent with radius Rmin , and the magnitude difference is 2.6 magnitude units. Even though there are 400 times as many
earthquakes of this size, they have little numerical significance.
Only ΔτS is varied when the goal is to find a model-matching statistical properties of a suite of strong-motion observations
from an earthquake. Increasing the value of ΔτS has two effects.
The first is that each subevent has a higher average slip, so theamplitude of the radiation from each subevent is increased.
Countering that trend, the number of subevents needed tomatch M 0 decreases. Because M 0 is not changed, the low-frequency radiation of the fault is not affected by ΔτS. However,the net effect of increasing ΔτS is an increase of the amplitude of
the high frequencies in the seismogram. Beyond the variability
in the locations and sizes of subevents, some source variability
arises when a random realization does not match the target seismic moment. In this case, the moment is matched by adjustingΔτS in equation (1) without changing the number, location, or
radii of subevents. Anderson (1997) found that ΔτS is proportional to radiated energy and apparent stress, so the CSM has theflexibility to match both of the important stress parameters.
Figure 1 illustrates the source in three ways: the point location of subevents with contours showing rupture time, subevent locations with circles of radius Ri showing the part of the
rupture represented by each subevent, and with contours illustrating the spatial distribution of slip represented by the sum of
all subevents. Figure 2 shows the moment rate represented bythis model and its displacement and acceleration spectra. In
Figure 2, the source time function is slightly longer than 11 s,
in which 11 s is the largest rupture time contoured in Figure 1.
The corner frequency corresponding to Rmin is about 3 Hz,
and the acceleration spectrum of this model is essentially flatat higher frequencies. The acceleration spectrum also shows
two corner frequencies, at about 0.2 Hz and about 3 Hz.
Green’s Function
The Green’s function is calculated using the approach of Lucoand Apsel (1983), using a FORTRAN code written by Y. Zeng.
Luco and Apsel derive the response of a layered half-space using
propagator matrices in the frequency–wavenumber domain,

15
0

10

20

30

40

0

10

20

30

40

10

20

30

40

(b) 0
5
10
15

(c) 0
5
10
15
0

▴ Figure 1. Composite source illustrated in three ways to emphasize different features. (a) Locations of the 528 subevents for one
realization of an M w 6.9 earthquake on a fault that is 40 km long and17.5 km wide. The dashed lines indicate the locations of the rupture
front in 1 s intervals. Radiation of each subevent is triggered when
the rupture front reaches its center. (b) The source radiation fromeach subevent is visualized as originating from a crack bounded bya circle centered on each subevent. The subevent radii are chosenrandomly from the distribution in equation (4). Subevents are notallowed to overlap the fault borders. (c) The slip modeled in eachsubevent is displayed here with an elliptical cross section for eachsubevent, superimposed and contoured.

organized in a way that has very good numerical stability.
The inverse transform over wavenumber is performed by anumerical integration. The time series is then recovered using
an inverse fast Fourier transform. Static offset is included in thesynthetic seismograms. It is possible to implement a frequencydependent Q in the Luco and Apsel (1983) model, but in thepresent version Q is independent of frequency. The treatmentof the frequency dependence observed in Q is described subsequently. A typical application calculates the synthetic Green’s
functions from 0 to 25 Hz. For the SCEC Broadband implementation, all seismograms are found for a common site condition,
but for other applications the shallow layers in the velocity model
can be tailored to the local site.

Seismological Research Letters−104−

Volume 86, Number 1

January/February 2015

69

10

–2
0.6

Q–1

Q(f) = 143 f
Approximation: 1/Q(f) = 1/1148 + 1/(120 f)

10

–3

0

5

10

15

20

f (Hz)▴ Figure 2. (a) Moment rate function of the source in Figure 1.
Units are in dyn·cm=s. (b) Fourier amplitude spectrum of moment
rate function. The Fourier amplitude of displacement seismograms
is similar but depends on azimuth because of radiation pattern anddirectivity. Units are in dyn·cm. Seismic moment of the source,
2:9 × 1026 dyn·cm, can be read from the intercept at f  0:01 Hz.
(c) Smoothed acceleration spectrum corresponding to the displacement spectrum in the center plot.
The Model for Q
The model for Q needs to satisfy several criteria. First, it should
be consistent with Q observations that find Q increasing withfrequency at high frequencies. Second, amplitudes of seismic
waves should have the same distance dependence as data.
Third, the S-wave Q and coda Q , measured from the synthetic
seismograms, should match observations. Finally, the trend of
the high-frequency shape should behave as exp−πκf  (Anderson and Hough, 1984). Considering that the spectrum of thesource model has a high-frequency asymptote of f −2 in displacement, or f 0 in acceleration (Fig. 2), a frequency-independent Q can cause this dependence. Also, I have modeled the
distance dependence of κ asdκκr  κ 0  κ~ r ≈ κ 0  rdr6(Anderson, 1986). The linear approximation in equation (6)was proposed by Anderson and Hough (1984) as a first approximation to the distance dependence.
Many regional Q models at high frequencies are modeledby a function of the form:
Q f   Q 0 f γ ;7

in which Q 0 is the value of Q at 1 Hz, and γ is usually inthe range 0 < γ < 1. For instance, in the SCEC Broadband
Platform simulations for California, the default Q structure
is given by Bayless:
70

Seismological Research Letters

Volume 86, Number 1▴ Figure 3. Example of the approximation to the Q model given byequation (9).
Q SCEC f   41  34βf 0:6 ;8

in which β is the shear velocity in the layer. For 0 ≤ γ ≤ 1,equation (7) can be approximated as
1
1
1
≈ 
Q 0f γ Q i Q d f9(e.g., Anderson, 1986). Q i and Q d in equation (9) can be foundso that the mismatch is under about 10% for the range of 2 ≤f ≤ 25 Hz (e.g., Fig. 3). Often, as in Figure 3, Q i is very large.
These high values are needed to model observations of dκ=dr(e.g., Hough and Anderson, 1988).
To understand the effect of this approximation, consider
the attenuation of a plane wave propagating a distance x for themodel of Q given in equation (9):


 


πxf
πxf
πx
≈ exp −exp −: 10
Ax  exp −γ
Q 0f βQ iβ
Q dβ
The first term on the right side causes the high frequencies todecay, and one can associate κ  x=Q i β. The second term on
the right side has no effect on the spectral shape, but it reduces
the amplitude uniformly across the spectrum.
To model κ0 , one needs low values of Q in the shallow
sediments. An initial assignment is given in Table 1. To matcha target value of κ0 , the values of Q for the upper 3 km areadjusted based on the slope of the Fourier SH response to
an impulse at 5 km depth. These low Q values in the shallow
crust rapidly damp any short-period surface waves so that evenat moderate distances the synthetic seismograms are dominatedby near-vertical propagating S waves.
Velocity Model Example
Figure 4 shows the velocity and Q model used for northern
California as modified from Goulet et al. (2015). All layers
January/February 2015

−105−

Table 1
Assignment of Q i for S Waves in Shallow Layers, Based on
Assimaki et al. (2008)
Shear Velocity for
Layers Shallower
Than 30 m (km=s)
0
0.1
0.5
1.0
2.0
10.0

Qi
5
5
7
10
10
10

Depth Rangez b ≤ 0:1 km
0:1 km ≤ z b ≤ 0:5 km
0:5 km ≤ z b ≤ 1:0 km
1:0 km ≤ z b ≤ 2:0 km
2:0 km ≤ z b ≤ 5:0 km

Qi
20
25
30
50
100

For the upper 30 m, Q i is interpolated from the first twocolumns. From 30 m to 3 km, Q i is determined from the depthz b of the bottom of the layer (third and fourth columns). Q i isgenerally defined in equation (9), but for these depths Q istreated as independent of frequency.

shallower than 30 m have been consolidated into a single layer
with 30 m thickness. Below the Moho, several 5 km thick layers
are introduced to bring the base of the model to 75 km depth;
velocity in each of these layers is increased using the Earth flattening approximation for SH waves (Aki and Richards, 2002),
βf z  a=a − zβz, in which a  6371 km is the radiusof the Earth and z is the layer depth. The P-wave velocity is
increased proportionately. The intent is to allow some of theenergy that is refracted into the mantle to return to the crustal
waveguide. Also note the high values of Q i in the crustal waveguide and low values in the upper 5 km.
One expects a trade-off between the velocity model andQ f . A result of addressing this trade-off is that Figure 4
0

0

5

5

has a small velocity gradient in the crustal waveguide createdby replacing the two thick original layers between 9 and 25 kmwith eight layers of 2 km thickness each, thus providing asmoother approximation to a velocity gradient in this depth
range. Without this gradient, the synthetic seismograms attenuate too rapidly. An alternative way to overcome that systematic effect is to increase Q more in those layers (to ∼4000),
but when that is done the distance dependence of κ does notmatch the trend of the measurements, at least for the Loma
Prieta earthquake, as discussed subsequently.
Coda
The coda in this implementation is calculated using the singlescattering model of Aki and Chouet (1975). It is applied as a
convolution with the Green’s functions. The coda function begins with a filtered impulse in which the filter incorporates theeffect of Q d from equation (10). It also adds reflections fromrandomly placed scatterers in the crustal waveguide. Figure 5
shows one of the Green’s functions for one station without andwith the coda. In this example, and in general, with the coda
the Green’s functions are beginning to take the realistic appearance of a small earthquake.

EXAMPLE: LOMA PRIETA CALCULATION
Synthetic Seismograms
Figures 6 and 7 show, respectively, the unfiltered synthetic acceleration, velocity, and displacement at station 0734, and thepseudorelative response spectra. The reader may judge theirrealism.
The average ground-motion amplitudes over multiple source
models are fairly stable, but individual realizations show variability. The standard deviation, say σ fault, is due to the variability of
0.4

10

10

15

15

20

20

25

25

30

30

Normalized Greens Function

Depth (km)0.2
Original

0
–0.2
–0.4
–0.6
–0.8
Coda Added

–1
–1.2
–1.4

35

35
Original
CSM

40

0

2

β

4

km/s

Qα
Qβ

α
6

8

40

0

500

–1.6
1000
Q

1500

0

20

▴ Figure 4. Original and modified velocity and Q model fornorthern California. Original models are shown with narrow lines,
and model used for the CSM are shown with heavy lines. α and βrefer to P and S waves, respectively.

40

60

80

t (s)

2000

▴ Figure 5. Example of a Green’s functions for station 0734 in theLoma Prieta earthquake. This is the north component for a fault
element near the south end of the rupture. (Top) As computed forthe velocity and Q model in Figure 4. (Bottom) The same component after convolution with a coda scattering function.

Seismological Research Letters−106−

Volume 86, Number 1

January/February 201571

sn0734(a)sn0734

3

10"E": max=35.97

2

20 cm/s2"N": max=69.06

Acceleration Response (cm/s2)

10

Z: max=51.65

1

10

0

10

−1

10

−2

10

−3

10

−2

−1

10

0

10

1

10
T (s)

2

10

10

20 s

sn0734(b)"E": max=6.097

10 cm/s"N": max=6.514

Z: max=5.387

20 s

sn0734(c)▴ Figure 7. East and north components of the pseudoacceleration response spectra (5% damping) for 50 realizations at station
0734. Black solid lines show averages for the two components,
and dashed lines show  one standard deviation.
the fault alone. For this Loma Prieta example, for 50 realizationsof peak acceleration at each station, σ fault decreases from ∼0:36 atshort distances to ∼0:24 at ∼80 km. For peak velocity, σ fault decreases from ∼0:32 to ∼0:2 in the same distance range. In the
future, controlling the variability of ΔτS can be used to control
the variability of σ fault .
Figure 8 compares data with the values of κ measured fromwhole-record synthetics using this model. Whole-record spectra are used because direct P waves overlap direct S waves atmany of the stations, so the whole-record window is a straightforward way to measure the same thing at all distances; the"E": max=6.227
0.1
Model
Data

0.09
0.08
0.07

3 cm

kappa (s)"N": max=7.107

0.06
0.05
0.04
0.03

Z: max=4.169

0.02
0.01
0

20 s

▴ Figure 6. Unfiltered synthetic realization of (a) acceleration,(b) velocity, and (c) displacement at station 0734. Notice the small
static offset in the horizontal displacements.

72

Seismological Research Letters

Volume 86, Number 1

0

10

20

30

40
50
r (km)

60

70

80

90

flt

▴ Figure 8. Whole-record measurements of κ for both observed
seismograms (open points, dashed line) and corresponding synthetic
seismograms (black points, solid line). The lines both have slope∼2:4 × 10−4 s= km. Observed spectra are not adjusted for V S30 .

January/February 2015−107−

data without violating the observed slope (dκ=dr). A full exploration of this trade-off is ongoing.
Calibration
When the average amplitudes of the CSM have a nonzero
residual for an event, the residual can be adjusted by changing
ΔτS . Figure 9 shows an example. The high-frequency residuals
decrease as the subevent stress drop increases. The lower frameplots the average residual for peak acceleration, taken to represent the high frequencies in general for the sake of this example, as a systematic function of ΔτS . Generating plots likethis for individual earthquakes identifies the optimal subevent
stress drop for each. With measurements of ΔτS from multiple
earthquakes, appropriate values can be determined for futureearthquakes.
At present, the SCEC calibration has not run enough cases
to determine if there is a dependence of ΔτS on magnitude,
focal mechanism, or region. At present, for California events
the CSM uses ΔτS  75 bars, even though ΔτS  100 bars
would be optimal for Loma Prieta.

SUMMARY▴ Figure 9. Loma Prieta summary results. Each of the top threeframes uses a different target subevent stress drop to find amplitudes at each station. The light lines show the corresponding residuals averaged over 50 realizations. The heavy lines show theaverage over the 40 stations. This figure shows how the residuals
for the response spectrum change as the subevent stress drop isincreased. The bottom frame shows how the residual for peak acceleration (T  0:01 s) depends on the value of ΔτS .
presence of the coda in these windows has increased κ above
the targeted value of κ 0  0:04 s for the S waves. The modelshows an increase of κ with distance similar to the data. Anderson (1991) averaged kappa estimates at several stations insouthern California, and an average trend through that modelfor κ~ r has the slope 2:9 × 10−4 s=km, which is similar to theslope in Figure 8. The distance dependence in Figure 8 is also
predicted by ray theory. As mentioned earlier, the slope(dκ=dr) in Figure 8 is controlled by Q i in the deeper parts ofthe crust. The attenuation is controlled by Q i, Q d , and thevelocity model. The introduction of the gradient in the lower
crust in Figure 4 achieved the right attenuation for Loma Prieta

This paper has summarized the theory behind the CSM forgenerating synthetic seismograms and given one example ofthe synthetics that are generated. For version 13.6 of the SCEC
Broadband Platform, used for 2013 gauntlet, the CSM was below the target range for GMPEs at 50 km. The new version,
with all of the features described in this article, had not beenformally evaluated as part of the version 14.3 gauntlet. Featuresdescribed here that were added subsequent to version 13.6 include the gradients in the crust and mantle, which are physically motivated adjustments to the velocity model, the additionof the coda, and the systematic evaluation of residuals and κ asa function of distance. Rules for the parameters controllingeach of these features will be incorporated into the next version
on the Broadband Platform.

ACKNOWLEDGMENTS
I thank D. Dreger and an anonymous reviewer for useful comments that improved this manuscript. This research was supportedby the Southern California Earthquake Center (SCEC). SCEC isfunded by NSF Cooperative Agreement EAR-1033462 and U.S.
Geological Survey Cooperative Agreement G12AC20038. The
SCEC contribution number for this paper is 1997.

REFERENCES
Aki, K., and B. Chouet (1975). Origin of coda waves: Source, attenuationand scattering effects, J. Geophys. Res. 80, 3322–3342.
Aki, K., and P. G. Richards (2002). Quantitative Seismology, Second Ed.,
University Science Books, Sausalito, California, 700 pp.
Anderson, J. G. (1986). Implication of attenuation for studies of theearthquake source, Earthquake Source Mechanics, Geophysical
Monograph 37, Maurice Ewing Series 6, American GeophysicalUnion, Washington, D.C., 311–318.

Seismological Research Letters−108−

Volume 86, Number 1

January/February 201573

Anderson, J. G. (1991). A preliminary descriptive model for the distancedependence of the spectral decay parameter in southern California,
Bull. Seismol. Soc. Am. 81, 2186–2193.
Anderson, J. G. (1997). Seismic energy and stress drop parameters for acomposite source model, Bull. Seismol. Soc. Am. 87, 85–96.
Anderson, J. G., and S. Hough (1984). A model for the shape of theFourier amplitude spectrum of acceleration at high frequencies,
Bull. Seismol. Soc. Am. 74, 1969–1994.
Anderson, J. G., and G. Yu (1996). Predictability of strong motions fromthe Northridge, California, earthquake, Bull. Seismol. Soc. Am. 86,
S100–S114.
Assimaki, D., W. Li, J. H. Steidl, and K. Tsuda (2008). Site amplification
and attenuation via downhole array seismogram inversion: Study ofthe 2003 Miyagi-Oki aftershock sequence, Bull. Seismol. Soc. Am.
98, 301–330.
Brune, J. N. (1970). Tectonic stress and spectra of seismic shear waves
from earthquakes, J. Geophys. Res. 75, 4997–5009.
Brune, J. N. (1971). Correction. Tectonic stress and spectra of seismic
shear waves from earthquakes, J. Geophys. Res. 76, 5002.
Frankel, A. (1991). High-frequency spectral falloff for earthquakes, fractal dimension of complex rupture, b-value, and the scaling ofstrength on faults, J. Geophys. Res. 96, 6291–6302.
Goulet, C., N. A. Abrahamson, P. G. Somerville, and K. E. Wooddell(2015). The SCEC Broadband Platform validation exercise forpseudo-spectral acceleration: Methodology for code validation in
the context of seismic hazard analyses, Seismol. Res. Lett. 86,
no. 1, doi: 10.1785/0220140104.
Hartzell, S., A. Frankel, P. Liu, Y. Zeng, and S. Rahman (2011). Model andparametric uncertainty in source-based kinematic models of earthquake ground motion, Bull. Seismol. Soc. Am. 101, 2431–2452.
Hough, S. E., and J. G. Anderson (1988). High-frequency spectra observed at Anza, California: Implications for Q structure, Bull. Seismol. Soc. Am. 78, 692–707.
Kanamori, H., and D. L. Anderson (1975). Theoretical basis of some empirical relations in seismology, Bull. Seismol. Soc. Am. 65, 1073–1095.
Khattri, K. N., G. Yu, J. G. Anderson, J. N. Brune, and Y. Zeng (1994).
Seismic hazard estimation using modelling of earthquake strong
ground motions: A brief analysis of 1991 Uttarkashi earthquake,
Himalaya and prognostication for a great earthquake in the region,
Curr. Sci. 67, 343–353.

74

Seismological Research Letters

Volume 86, Number 1

Luco, J. E., and R. J. Apsel (1983). On the Green’s function for a layered
half-space, part I, Bull. Seismol. Soc. Am. 73, 909–929.
Sato, T., and T. Hirasawa (1973). Body wave spectra from propagating
shear cracks, J. Phys. Earth 21, 415–431.
Su, F., J. G. Anderson, S.-D. Ni, and Y. Zeng (1998). Effect of site amplification and basin response on strong motion in Las Vegas, Nevada, Earthq. Spectra 14, 357–376.
Su, F., J. G. Anderson, and Y. Zeng (1998). Study of weak and strong
ground motion including Nonlinearity from the Northridge, California, earthquake sequence, Bull. Seismol. Soc. Am. 88, 1411–1425.
Yu, G. (1994). Some aspects of earthquake seismology: Slip partitioning
along major convergent plate boundaries; composite source modelfor estimation of strong motion; and nonlinear soil response modeling, Ph.D. Thesis, University of Nevada, Reno, 144 pp.
Yu, G., K. N. Khattri, J. G. Anderson, J. N. Brune, and Y. Zeng (1995).
Strong ground motion from the Uttarkashi, Himalaya, India earthquake: Comparison of observations with synthetics using thecomposite source model, Bull. Seismol. Soc. Am. 85, 31–50.
Zeng, Y., and J. G. Anderson (1996). A composite source model of the1994 Northridge earthquake using genetic algorithms, Bull. Seismol.
Soc. Am. 86, no. 1B, S71–S83.
Zeng, Y., and C. Chen (2001). Fault rupture process of the 20 September
1999 Chi-Chi, Taiwan, earthquake, Bull. Seismol. Soc. Am. 91,
1088–1098.
Zeng, Y., J. G. Anderson, and F. Su (1995). Subevent rake and random
scattering effects in realistic strong ground motion simulation, Geophys. Res. Letters 22, 17–20.
Zeng, Y., J. G. Anderson, and G. Yu (1994). A composite source model
for computing realistic synthetic strong ground motions, Geophys.
Res. Lett. 21, 725–728.
Zeng, Y., F. Su, and K. Aki (1991). Scattered wave energy propagation in arandom isotropic scattering medium, J. Geophys. Res. 96, 607–619.

John G. Anderson
University of Nevada
Reno, Nevada 89557 U.S.A.
jga@unr.edu

Published Online 17 December 2014

January/February 2015−109−

Bulletin of the Seismological Society of America, Vol. 107, No. 6, pp. 2561–2577, December 2017, doi: 10.1785/0120160361Ⓔ

Fault-Scaling Relationships Depend on the Average Fault-Slip Rate
by John G. Anderson, Glenn P. Biasi,* and Steven G. Wesnousky

Abstract

This study addresses whether knowing the slip rate on a fault improves
estimates of magnitude (M w ) of shallow continental surface-rupturing earthquakes.
Based on 43 earthquakes from the database of Wells and Coppersmith (1994), Anderson et al. (1996) suggested previously that the estimates of M w from rupture length(LE ) are improved by incorporating the slip rate of the fault (SF ). We re-evaluate thisrelationship with an expanded database of 80 events, which includes 56 strike-slip, 13
reverse-, and 11 normal-faulting events. When the data are subdivided by fault mechanism, magnitude predictions from rupture length are improved for strike-slip faults
when slip rate is included but not for reverse or normal faults. Whether or not theslip-rate term is present, a linear model with M w ∼ log LE over all rupture lengths
implies that the stress drop depends on rupture length—an observation that is notsupported by teleseismic observations. We consider two other models, including one
we prefer because it has constant stress drop over the entire range of LE for any constant value of SF and fits the data as well as the linear model. The dependence on slip
rate for strike-slip faults is a persistent feature of all considered models. The observed
dependence on SF supports the conclusion that for strike-slip faults of a given length,
the static stress drop, on average, tends to decrease as the fault-slip rate increases.

Electronic Supplement: Table of earthquakes and parameters.

Introduction
Models for estimating the possible magnitude of anearthquake from geological observations of the fault length
are an essential component of any state-of-the-art seismichazard analysis. The input to either a probabilistic or deterministic seismic-hazard analysis requires geological constraints because the duration of instrumental observations ofseismicity is too short to observe the size and to estimate theoccurrence rates of the largest earthquakes (e.g., Allen, 1975;
Wesnousky et al., 1983). Thus, wherever evidence in thegeological record suggests earthquake activity, it is essential
for the seismic-hazard analysis to consider the hazard fromthat fault, and an estimate of the magnitude of the earthquake(Mw ) that might occur on the fault is an essential part of theprocess. The primary goal of this study is to determine if
magnitude estimates that are commonly estimated from fault
length (LE ) can be improved by incorporating the slip rate(SF ) of the fault.
Numerous models for estimating magnitude from rupture
length have been published. Early studies were by Tocher(1958) and Iida (1959). Wells and Coppersmith (1994) published an extensive scaling study based on 244 earthquakes.
*Now at U.S. Geological Survey, 525 S. Wilson Avenue, Pasadena,
California 91106; gbiasi@usgs.gov.

Some of the more recent studies include Anderson et al.
(1996), Hanks and Bakun (2002, 2008), Shaw and Wesnousky(2008), Blaser et al. (2010), Leonard (2010, 2012, 2014), andStrasser et al. (2010). For probabilistic studies and for earthquake source physics, it is valuable to try to reduce the uncertainty in these relations. Motivated by Kanamori and Allen(1986) and Scholz et al. (1986), Anderson et al. (1996; hereafter, AWS96) investigated whether including the fault-slip
rate on a fault improves magnitude estimates given rupture
length. They found that it does, and proposed the relationshipMw  5:12  1:16 log LE − 0:20 log SF , thus indicating thatslip rate is a factor. A physical interpretation of a significant
dependence on slip rate is that, for a common rupture length,faults with higher slip rates tend to have smaller static stressdrop. Since the publication of AWS96, the number of earthquakes with available magnitude, rupture length, and slip-rate
estimates has approximately doubled. This article considers
whether these new data improve or modify the conclusions
from the earlier study.
One consideration in developing a scaling model is thatseismological observations have found stress drop in earthquakes to be practically independent of magnitude. Kanamori and Anderson (1975) is one of the early papers tomake this observation. Recent studies that have supported

2561

−111−

2562

J. G. Anderson, G. P. Biasi, and S. G. Wesnousky

Table 1
Earthquakes from 1968–2011 Used in This Study
Event
NumberEvent Name

Event Date(yyyy/mm/dd)

Mw

RuptureLength (km)Slip Rate(mm/yr)Mechanism*

2
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
25
26
28
29
30
31
32
33
34
35
36
37
38
40
41
42
43
44
45
46
47

Fukushima-Hamadori, Japan
Yushu, China
El Mayor–Cucapah
Wenchuan, China
Kashmir, Pakistan
Chuya, Russia (Gorny Altai)
Denali, Alaska
Kunlun, China
Düzce, Turkey
Hector Mine, California
Chi-Chi, Taiwanİzmit, Turkey
Fandoqa, Iran
Manyi, China
Sakhalin Island (Neftegorsk), Russia
Northridge, California
Landers, California
Luzon, Philippines
Rudbar, Iran
Loma Prieta, California
Superstition Hill, California
Edgecumbe, New Zealand
Marryat, Australia
Morgan Hill, California
Borah Peak, Idaho
Coalinga, California
Sirch, IranCorinth, Greece
Corinth, Greece
Daofu, China
El Asnam (Ech Cheliff), Algeria
Imperial Valley, California
Coyote Lake, California
Tabas, Iran
Bob-Tangol, Iran
Motagua, Nicaragua
Luhuo, China
San Fernando, California
Tonghai, China
Dasht-e-Bayaz, Iran
Borrego Mtn, California

2011/04/11
2010/04/14
2010/04/04
2008/05/12
2005/10/08
2003/09/27
2002/11/03
2001/11/14
1999/11/12
1999/10/16
1999/09/21
1999/08/17
1998/03/14
1997/11/08
1995/05/27
1994/01/17
1992/06/28
1990/07/16
1990/06/20
1989/10/17
1987/11/24
1987/03/02
1986/03/03
1984/04/24
1983/10/28
1983/05/02
1981/07/29
1981/02/25
1981/03/04
1981/01/24
1980/10/10
1979/10/15
1979/08/06
1978/09/16
1977/12/19
1976/02/04
1973/02/06
1971/02/09
1970/01/04
1968/08/31
1968/04/09

6.7
6.8
7.3
7.9
7.6
7.2
7.8
7.7
7.1
7.1
7.7
7.5
6.6
7.4
7.0
6.7
7.2
7.7
7.4
6.8
6.6
6.4
5.8
6.1
6.9
6.4
7.1
6.1
5.9
6.7
6.9
6.4
5.8
7.4
5.8
7.5
7.5
6.8
7.2
7.1
6.6

15
52
117
240
70
70
340
450
40
48
72
145
22
170
40
21
77
112
80
35
25
15.5
13
20
36
25
65
14
15
44
36
36
14
90
19.5
230
90
19
60
80
33

0.02
12
2.5
1.3
3.1
0.5
12.4
10
15
0.6
12.9
12
2
3
4
0.4
0.4
15
1
3.2
3
2
0.005
5.2
0.15
1.4
4.3
1.7
0.3
12
0.8
17
11.9
1.3
4
12
14
1.8
2
5
6.7

N
S
S
R
R
S
S
S
S
S
R
S
S
S
S
R
S
S
S
R
S
N
R
S
N
R
S
N
N
S
R
S
S
R
S
S
S
R
S
S
S

*N, normal, S, strike slip, R, reverse.

this result include Allmann and Shearer (2009) and Baltay
et al. (2011). Apparent exceptions have been reported basedon Fourier spectra of smaller earthquakes, but as magnitudedecreases, attenuation can cause spectral shapes to behave
the same as they would for decreasing stress drop (e.g., Anderson, 1986). Studies that have taken considerable care to
separate these effects have generally concluded that the average stress drop remains independent of magnitude down to
extremely small magnitudes (e.g., Abercrombie, 1995; Ideet al., 2003; Baltay et al., 2010, 2011). However, all of thesestudies find that for any given fault dimension, the range ofmagnitudes can vary considerably (e.g., Kanamori and Allen,1986). Despite this variability, it seems reasonable to evaluate a scaling relationship that is based on a constant stressdrop before considering the additional effect of the fault-slip
rate. This vision guides the development of the consideredscaling relationships. Details of these models for the relationship of stress drop and the fault dimensions are deferred to
the Appendix. The following sections describe the data,
present the summary equations for three alternative models,
fit the alternative models to the data, and discuss the results.

Data
Tables 1 and 2 give the preferred estimates of M w , LE ,and SF for the earthquakes used in this analysis. These values
and their corresponding uncertainty ranges are given in
Table S1 (Ⓔ available in the electronic supplement to this−112−

Fault-Scaling Relationships Depend on the Average Fault-Slip Rate

2563

Table 2
Earthquakes from 1848–1967 Used in This Study
Event
NumberEvent Name

Event Date(yyyy/mm/dd)

Mw

RuptureLength (km)Slip Rate(mm/yr)Mechanism*

48
49
51
52
53
54
55
56
57
58
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88

Mudurnu Valley, Turkey
Parkfield, California
Alake Lake or Tuosuohu Lake or Dulan, China
Ipak or Buyin-Zara, Iran
Hebgen Lake, Montana
Gobi-Altai, Mongolia
San Miguel, Mexico
Fairview Peak, Nevada
Dixie Valley, Nevada
Yenice–Gonen, Turkey
Gerede-Bolu, Turkey
Tosya, Turkey
Tottori, Japan
Niksar-Erbaa, Turkey
Imperial Valley, California
Erzincan, Turkey
Tuosuo Lake, Huashixia, China
Parkfield, California
Long Beach, California
Changma, China
Fuyun, China
North Izu, Japan
Laikipia, Kenya
Tango, Japan
Luoho-Qiajiao (Daofu), China
Haiyuan, China
Pleasant Valley, Nevada
Chon-Kemin (Kebin), Kazakhstan
San Francisco, California
Bulnay, Mongolia
Laguna Salada, Mexico
Rikuu, Japan
Nobi/Mino-Owari, Japan
Canterbury, New Zealand
Sonora, Mexico
Owens Valley, California
Hayward, California
Fort Tejon, California
Marlborough, New Zealand

1967/07/22
1966/06/28
1963/04/19
1962/09/01
1959/08/18
1957/12/04
1956/02/14
1954/12/16
1954/12/16
1953/03/18
1944/02/01
1943/11/26
1943/09/10
1942/12/20
1940/05/19
1939/12/25
1937/01/07
1934/06/08
1933/03/10
1932/12/25
1931/08/10
1930/11/25
1928/01/06
1927/03/07
1923/03/24
1920/12/16
1915/10/03
1911/01/03
1906/04/18
1905/07/23
1892/02/23
1896/08/31
1908/06/28
1888/09/01
1887/05/13
1872/03/26
1868/10/01
1857/01/09
1848/10/16

7.3
6.2
7.0
7.0
7.3
8.1
6.6
7.1
6.6
7.3
7.3
7.6
6.9
6.8
7.1
7.8
7.6
6.2
6.4
7.6
7.9
6.9
6.8
7.0
7.3
8.0
7.3
8.0
7.9
8.5
7.2
7.2
7.4
7.1
7.2
7.4
6.9
7.8
7.5

80
28
40
100
25
260
20
46
47
60
155
275
33
50
60
330
150
25
22
149
160
28
38
35
80
237
61
177
497
375
42
40
80
65
101.8
110
61
339
134

18
30
12
1
0.5
1
0.3
0.14
0.5
6.8
18
19
0.3
19
17
19
11
30
1.1
5
0.3
2.4
0.18
0.3
10
7
0.1
2
21
3
2.5
1
1.6
14
0.08
3.5
8
25
5.6

S
S
S
R
N
S
S
N
N
S
S
S
S
S
S
S
S
S
S
S
S
S
N
S
S
S
N
R
S
S
S
R
S
S
N
S
S
S
S

*N, normal, S, strike slip, R, reverse.

article), along with the references for all estimates. Events
considered for analysis come from AWS96 and Biasi andWesnousky (2016). Some AWS96 events were not used because uncertainties in one or more of the magnitude, length,
or slip-rate parameters were considered too large or too
poorly known to contribute to the parametric regressions.
Events in Biasi and Wesnousky (2016) were selected on
the basis of having a well-mapped surface rupture and nongeologically estimated magnitude. Their list builds on the listof fault ruptures of Wesnousky (2008) by adding more recent
events and by including surface ruptures newly documented
by geologic field work. Interested readers are referred tothese previous papers for further description of each event.
Overall, the database is heavily weighted toward surfacerupturing earthquakes. Some events in Biasi and Wesnousky(2016) were not included for lack of a resolved fault-slip rate,

or because their rupture lengths were too short. Events with
LE < 15 km were generally not included. The smallest preferred estimate of M w is 5.7.
Earthquakes after 1900 were only included if some independent (nongeologic) means was available to estimate
magnitude. Moment estimates from waveform modelingwere preferred to body-wave magnitudes where both wereavailable. The six earthquakes prior to 1900 are particularly
well documented, as described in Ⓔ the electronic supplement. Because LE is known for these events and the uncertainty in M w introduced by uncertain depth of faulting is lessthan 0.1, the measured slip in these events controls M w . It
follows that estimating M w from LE alone for these events
is not circular. The rupture length is normally taken as the
distance between the ends of primary coseismic surface rupture. The sum of the lengths of overlapping traces may be

−113−

2564

J. G. Anderson, G. P. Biasi, and S. G. Wesnousky

80

10

0

All
Strike slip
ReverseNormal70

Exceedance rate, yr –1

Cumulative Number

60
50
40
30

10

–1

10

–2

20
10
0
1800

1850

1900

1950

2000

2050

5.5

Year

6

6.5

7

7.5

8

8.5

9

Mw

Figure 1.

Figure 2.

used as the length in the analysis (e.g., event 53, Hegben
Lake, 1959) where the overlapping portions were judged tocontribute materially to the moment release. Rupture lengths
based on aftershock distributions have generally been
avoided, with the exception of six moderate strike-slip
events, all in California. These were retained for continuity
with AWS96 and for support of the regressions at moderatemagnitudes. None control the results. Fault-slip rates aretaken from offsets of geologic features 10–100 ka in age,where possible, to represent a stable recent slip-rate estimate.
Fault-slip rates from paleoseismic offsets of one or a few
individual earthquakes were avoided, because it is not clear
how that activity would relate to the longer term average slip
rate. Similarly, fault-slip rates from geodetic estimates were
avoided where possible because they measure the currentday rate but may not represent the longer term average. Fault
creep effects were considered, but no corrections were attempted in the database. First, creep is believed to affect only
a few percent or less of events, and at a fraction of the full slip
rate. Second, uncertainty in fault-slip rate will be seen below
to have little effect on the regression.
Figure 1 shows the cumulative number of earthquakes
used as a function of time. From 1954 to 2013, the rate ofusable events is relatively steady, about 0.9 events per year.
The rate is lower prior to ∼1954 suggesting that the earlier
historical record is less complete.
The earthquakes are separated into general categories ofstrike-slip, normal, and reverse faulting. Figure 2 shows theexceedance rates of considered earthquakes in each of thesecategories as a function of magnitude, both combined andseparated by focal mechanism. To estimate the rates, thenumber of earthquakes for each of the curves was divided
by 100 yrs. This is obviously an approximation, but considering Figure 1, the events prior to ∼1910 may roughly compensate for the missing events since 1910. For instance, Figure 2
suggests that continental events that cause surface rupture with

Mw ≥ 7:0 have occurred at a rate of about 0:5 yr−1 , or roughly
once every two years. The rates of strike-slip, reverse, andnormal mechanisms are about 0.4, 0.075, and 0:045 yr−1 .
Rounded to the nearest 5%, this implies that about 75% of
those events were strike slip, about 15% had reverse mechanisms, and about 10% had normal mechanisms.
Figure 3 shows maps with locations of all events, using
different symbols to distinguish among mechanisms. The insets show more details on locations of events from the western
United States, the eastern Mediterranean region, and Japan.
Figure 4 plots the preferred slip rates versus the preferred rupture lengths. Figure 4a emphasizes the overall distribution of the data, while 4b highlights the 56 strike-slip
faults, 4c highlights the 13 reverse faults, and 4d highlights
the 11 normal faults. The combined data in Figure 4a are
distinctly upper triangular. The points along the diagonal associate an increase in fault length with an increase in faultslip rate, which in turn is likely a function of cumulative slip
(e.g., Wesnousky, 1988, 1999) that does not depend on themechanism. The data above the diagonal show that (1) theentirety of long faults and fault systems does not always
break and that (2) small fast faults may exist. There are twoalternatives to explain the lack of long ruptures on faults withlow slip rates. The first could be purely statistical because
events in the lower right corner of the plot could be too rare
to be represented in the historical record. Alternatively, due
to what Perrin et al. (2016) call “the competition between
damage and healing processes,” faults with slow slip rates
might, during the interseismic period, be sufficiently affectedby differential healing, influences from adjacent faults, orother processes that long ruptures on slow faults never occur.
Figure 4b–d emphasizes the data available to search forslip-rate dependency for the three fault types. Figure 4b

Cumulative number of events used in this analysis(Tables 1 and 2), shown as a function of time. The color versionof this figure is available only in the electronic edition.

Event rates, as a function of magnitude and event
types. The rates are estimated based on the approximation thatthe data represent about 100 yrs of seismicity, as discussed in
the Data section. The color version of this figure is available onlyin the electronic edition.

−114−

Fault-Scaling Relationships Depend on the Average Fault-Slip Rate

Western United States

Northeastern Mediterranean

2565

Japan

Figure 3.

Locations of events considered in this study. Open circles show locations of events with a strike-slip mechanism. Triangles
represent reverse events, and inverted triangles represent normal mechanisms.

shows a distribution of points spanning rupture lengths
mostly between 20 and 400 km and slip rates between 0.3
and 30 mm=yr. Rupture lengths for reverse faults (Fig. 4c)range from 13 to 240 km although most are between 20 and
100 km. Slip rates for reverse faults are mostly between 0.4
and 4 mm=yr, with outliers at 0.005 and 12:9 mm=yr. Rupture lengths and slip rates for normal faults (Fig. 4d) range
from 14 to 102 km and 0.08 to 2 mm=yr, but are unevenly
distributed within these limits.

Modeling Approaches
The effect of slip rate is tested against three model
shapes for the scaling relationship to confirm that it is not
an artifact of a particular assumption for how magnitude depends on rupture length. The first M1 explores a linear regression of Mw with the logs of length and slip rate:
EQ-TARGET;temp:intralink-;df1;55;107

M w  c0  c1 log LE  c2 log

SF
;
S0

1

in which LE is the rupture length (measured along strike) of a
specific earthquake; M w is the reported moment magnitude
for the respective earthquake (Kanamori, 1977); SF is the slip
rate of the fault on which the earthquake occurred determined from geological observation; S0 is the average of thelogs of all slip rates in the data set being considered (e.g.,
strike-slip faults, normal faults, etc.); and c0 , c1 , and c2
are coefficients of regression to be determined. Mathematically, c0 trades off with −c2 log S0 , which allows the parameter S0 to be rounded to two significant digits. In this model,
setting SF  S0 is mathematically equivalent to settingc2  0 and thus also equivalent to the model approach usedby Wells and Coppersmith (1994) and others who estimate alinear dependence of M w on log LE without including theslip rate on the fault. Two misfit parameters are considered.
The first σ 1L is the standard deviation of the difference between observed and predicted magnitudes when c2  0, so
only LE is used to estimate Mw , while σ 1S is the corresponding standard deviation when the slip-rate term in equation (1)is incorporated. A consequence of the assumed model M1 is

−115−

2566

(a) 10 2

(b)

10 2

S F , mm/yr

S , mm/yr

J. G. Anderson, G. P. Biasi, and S. G. Wesnousky

10 0

F

10 0

10 –2

10 1

Strike slip
ReverseNormal

10 –2

10 2

10 3

10 1

L , km10 2

10 3

L , km(c) 10 2

(d)

10 2

S , mm/yr

E

S F , mm/yr

E

10 0

F

10 0

Strike slip
ReverseNormal

10 –2

10 1

10 2

Strike slip
ReverseNormal

10 –2

10 3

10 1

L , km10 2

10 3

L , kmE

E

Figure 4.

Rupture length–slip-rate distribution of the data in Tables 1 and 2. All points are shown combined in each part of the figure,(a) each fault type is given equal emphasis, (b) strike-slip faults are emphasized, (c) reverse faults are emphasized, and (d) normal faults areemphasized. The color version of this figure is available only in the electronic edition.

that unless c1 fortuitously equals 2/3, stress drop increases
for large earthquakes as a function of rupture length LE , regardless of whether slip rate is included or not (Table A1).
The second M2 constrains the slope to give constantstress drop for small and large earthquakes with a slopechange at the break-point magnitude M bp . The stress drop
for small and large earthquakes is allowed to differ:
 
 
8
< M bp  c1C log LE  c2 log SF LE < Lbp
Lbp  S0 
Mw 
2
: M bp  c1L log LE  c2 log SF LE ≥ Lbp ;Lbp
S0
EQ-TARGET;temp:intralink-;df2;55;200

in which c1C  2 and c1L  2=3 for rupture lengths that areless than or greater than Lbp , respectively, the rupture length
where slope changes from 2 to 2/3. The three unknown
parameters in model M2 are Lbp the rupture length where theslope changes from 2 to 2/3, M bp the magnitude at that transition, and c2 which is again the sensitivity of magnitude tofault-slip rate. Ruptures of length less than Lbp are consid-

ered to be a small earthquake and scale like a circular rupture
in Table A1, implying that constant stress drop occurs whenc1C  2. An earthquake with rupture length greater than Lbp
is considered to be a large earthquake and corresponds to oneof the models for a long fault in Table A1 (depending on fault
mechanism), for which the value c1L  2=3 results in constant stress drop. However, equation (2) does not require thestress drop for the small earthquakes to be the same as thestress drop for large earthquakes. Equation (2) has the same
number of unknown parameters to be determined from thedata as equation (1). The two standard deviations of the misfit for model M2 are σ 2L and σ 2S , which correspond directly
to the parameters σ 1L and σ 1S of model M1.
The third model M3 is derived from the model of Chinnery (1964) for a vertical strike-slip fault that ruptures thesurface. It is assumed that stress drop for the top center ofthe fault in this model ΔτC is constant across all rupture
lengths and magnitudes:−116−

Fault-Scaling Relationships Depend on the Average Fault-Slip Rate

EQ-TARGET;temp:intralink-;df3;55;733

2567



 
8
>
S
>
< 2 log LE  23 log ΔτC  23 log C2 2πCγ − 16:1  c2 log SF0
LW


 
Mw 
2
>
SF
max
>
: 23 log LE  23 log ΔτC  23 log 2πWCγ − 16:1  c2 log S0in which
EQ-TARGET;temp:intralink-;df5;313;664

EQ-TARGET;temp:intralink-;df4;55;652

cos γ sin γ3  4 sin γ
Cγ  2 cos γ  3 tan γ −:1  sin γ24

Details on the development of the model M3 equations areprovided in the Appendix. The value γ is the angle from the
top center of the fault to either of the bottom corners, that is,
tan γ  2W E =LE , in which W E is the down-dip width of theearthquake rupture. The model variables include four parameters. These are the aspect ratio of the fault for small ruptures
CLW  LE =W E , the stress drop ΔτC , the coefficient thatquantifies the slip-rate dependence c2 , and the maximum
fault width W max . Equation (3) assumes that the aspect ratio
is constant for small earthquakes and that when the selected
aspect ratio in combination with LE implies a width greaterthan W max the width is set to W max . For model M3, the twostandard deviations of the misfit are σ 3L and σ 3S , corresponding to the parameters σ 1L and σ 1S of model M1. As written,
the coefficients of the term in log LE appear to be the same as
in model M2, but for the long ruptures, γ depends on LE , so
the term with Cγ modifies the slope.
Model equations (1)–(3) require different strategies toobtain their unknown coefficients. The simplest way to findthe unknown coefficients for equation (1) is using a linear
least-squares regression, which minimizes the misfit of theprediction of M w but does not account for uncertainties in
LE or SF. AWS96 approached this difficulty by carryingout multiple regressions for points chosen at random within
the range of allowed values of all three parameters, and thenlooked at the distribution of derived values of the coefficientsof the regression. Alternative approaches to find the coefficients, described variously as “total least squares” or “general orthogonal regression” (e.g., Castellaro et al., 2006;
Castellaro and Bormann, 2007; Wikipedia article “Total
Least Squares,” see Data and Resources) were also considered for this analysis. The approach by AWS96 turned out togive the least biased results for a set of synthetic data with anuncertainty model that we considered to be realistic and consistent with the actual data, so their approach is also used inthis study. The parameters for equation (1) were determinedfrom 10,000 realizations of the randomized earthquake
parameters to find the distributions of coefficients.
In implementing the AWS96 approach, M w , LE , and SF
are chosen at random from the range of uncertainties given inⒺ the electronic supplement. The probability distributions forthe randomized parameters reflect that uncertainty ranges are
not symmetrical around the preferred value. The preferredvalue is set to be the median. As an example, the probability
distribution for the ith randomized value of LE is as follows:pLE  LE
CLW < W max
LE
CLW ≥ W max ;1minLpref
E −LE 1max
LE −Lpref
E caseAcaseB;

3

5

in which case A has probability of 0.5, case B has probability
max are the minimum and maximum of theof 0.5, Lmin
E and LErange on the rupture length, respectively, and Lpref
E is the preferred value. The seismic moment and slip rate are randomized
using the same algorithm, and M w is found from the randomized moment. The standard deviations of the misfit σ 1L
and σ 1S are the average values from the multiple realizations.
Equation 2 has the additional complication of being
nonlinear in Lbp . We approach the solution by reorganizing
equation (2) as
 
 
S
L
Mbp  c2 log F  M w − c1x log E ;
6
S0
Lbp
EQ-TARGET;temp:intralink-;df6;313;499

in which c1x is either c1C or c1L depending on LE . Assuming
a value for Lbp, it is straightforward to find the unknown coefficients Mbp and c2 . We considered a set of closely spaced
values of Lbp from the smallest to the longest rupture length
in the data, and choose the value with the smallest total misfit. For each trial value of Lbp , we solved for the unknown
coefficients 10,000 times with values of Mw , LE , and SF
randomized as in equation (5), and our preferred model isthe mean of the coefficients from the multiple realizations.
Model M3 (equation 3) has four unknown parameters, inwhich the effects of CLW and W max are nonlinear (Fig. A1).
For this reason, a grid of values of CLW and W max was
searched; there were 506 points on this grid. For each grid
point, ΔτC and c2 were determined by linear least squares for
10,000 randomly chosen realizations of Mw , LE , and SF . The
average values of ΔτC and c2 were found from the distributions of these realizations, together with the average values of
σ 3L and σ 3S . This permitted creating a contour plots of σ 3L
and σ 3S as a function of the trial values of CLW and W max .
The minima in σ 3L and σ 3S did not generally occur for thesame combinations of CLW and W max . Because the results ofmodel M3 might potentially be used for faults where slip rateis unknown, we minimized σ 3L . The minimum in σ 3L is
broad compared with the grid spacing of CLW and W max , so
the values that are used come as near as possible, within theminimum of σ 3L , to minimize σ 3S as well. The grid limits
considered the maximum fault widths from 10 to 20 km forstrike-slip faults, whereas for reverse and normal faulting thegrid limits considered the maximum fault widths from 18 to
30 km. The larger widths were considered because of thesuggestions of King and Wesnousky (2007), Hillers andWesnousky (2008), and Jiang and Lapusta (2016) that a dynamic rupture in a large earthquake might reasonably extend

−117−

2568

(a)

J. G. Anderson, G. P. Biasi, and S. G. Wesnousky(a) 9 Strike slip

Strike slip

Mw

Mw

8

7
SF > 4.8 mm/yr
SF ≤ 4.8 mm/yr

6
w

5
1
10

2L

L bp = 73.8; M bp = 7.38

= 0.238
10

2

10

3

10

2

10

3

10

2

10

3

LE , kmResidual from Top Plot
M w = –0.176 log (S F / 4.8 mm/yr)

1
w

Mw

Mw

0.5
0
–0.5
2S

–1
–1
10

10

0

10

= 0.215

1

S , mm/yr
F(b) 9 Reverse(b) 9 Reverse

8

Mw

Mw

8

7

7
SF > 1.1 mm/yr

SF > 1.1 mm/yr
SF

6
1L
5
1
10

M w = 5.12 + 1.15 log LE= 0.322
10

2

10

5
1
10

3

L bp = 46.4; M bp = 7.23

= 0.281
10

Residual from Top Plot
M w = 0.264 log ( SRF / 1.1 mm/yr )Residual from Top Plot
M w = 0.169 log ( SF / 1.1 mm/yr)

Mw

0.5

0

0
–0.5

–0.5
–1 1S
–3
10

= 0.238
10

–2

10

–1

10

0

10

1

10

2S

–1
–3
10

2

10

–2

10

F

(c) 9 Normal

010= 0.253

1

8

7

Mw

Mw

10

(c) 9 Normal

8

SF > 0.25 mm/yr
SF

6SF > 0.25 mm/yr

10 2

5
1
10

10 3

= 0.289
2L
10

2

Residual from Top Plot
M w = –0.107 log (SF / 0.25 mm/yr)

1

0.5

0.5

Mw

0
–0.5
1S
–1
–2
10

L bp = 24.3; M bp = 6.8

LE , kmResidual from Top Plot
M w = -0.115 log ( SRF / 0.25 mm/yr )1SF ≤ 0.25 mm/yr

6

M w = 5.25 + 1.02 log LE= 0.318
1L

7

0.25 mm/yr

LE , km

Mw–1

SF , mm/yr

S , mm/yr

5
10 1

2

1

1

Mw

2L

LE , km

LE , km0.5

SF ≤ 1.1 mm/yr

6

1.1 mm/yr

0–0.5= 0.303
10

–1

10

0

10

–1
10 -2

1

S , mm/yr

10 0= 0.277
10 1

SF , mm/yr

F

Figure 5.

2S

10 –1

Model M1 (equation 1) for (a) strike-slip, (b) reverse,and (c) normal faults. For each mechanism, the upper frame shows
Mw plotted as a function of LE . Points are all the preferred values, as
given in Tables 1 and 2. Solid points represent low slip-rate faults.
The solid line uses coefficients given in Table 3 for SF  S0. The
lower frame shows the residuals δMw of the points in the upper
frame from the solid line. The line in the lower frame shows thepredicted effect of SF based on the coefficients in Table 3, thatis, δM w  c2 logSF =S0 . For strike-slip faults, the significant effect of fault-slip rate is seen in the clear separation of low and highslip-rate faults in the upper panel, and the negative slope of the fit to
the residuals in the lower panel. For reverse and normal faults, thesparse data suggest a different trend in the residuals, indicating thatmixing the three mechanism types is not appropriate.

Figure 6.

Model M2 (equation 2) for (a) strike-slip, (b) reverse,and (c) normal faults. Coefficients for the lines are given in Table 4.
Other figure details are the same as in Figure 5.

deeper than the brittle crustal depths associated with microearthquakes.

Analysis Results
Figures 5–7 show results for models M1–M3, respectively, for strike-slip, reverse-, and normal-faulting earthquakes. For each mechanism, the curve in the upper frame^ w for SF  S0. The
shows predicted values of magnitude M−118−

Fault-Scaling Relationships Depend on the Average Fault-Slip Rate

2569(a) 9 Strike slip

Table 3
Coefficients for Model M1, Use in Equation (1), the Different
Fault Types Considered Separately, and Earthquakes Listed inTables 1 and 2

Mw

8

7
SF > 4.8 mm/yr
SF ≤ 4.8 mm/yr

6
3L
5
1
10

MC (

= 0.236

10

2

10

3

LE , kmResidual from Top Plot
Mw = –0.17 log(S F / 4.8 mm/yr )

1
0.5

Mw

Parameter

Strike Slip

ReverseNormalc0c1
c2
S0 (mm/yr)
σ 1L
σ 1S

4.73 ± 0.062
1.30 ± 0.031−0.198 ± 0.023
4.8
0.241
0.211

5.12 ± 0.11
1.15 ± 0.065
0.264 ± 0.036
1.1
0.322
0.238

5.25 ± 0.18
1.02 ± 0.12
−0.115 ± 0.109
0.25
0.318
0.303

=24.9 bars, CLD =3.8, W max =15) [SF =4.8 mm/yr]

0
–0.5
–1 3S
10 –1

= 0.214
10 0

10 1

10 2

SF , mm/yr

Table 4(b) 9 Reverse

Coefficients for Model M2, Use in Equation (2), the Different
Fault Types Considered Separately, and Earthquakes Listed inTables 1 and 2

Mw

8

7

Parameter

Strike Slip

ReverseNormalLbp
Mbp
c2
S0 (mm=yr)
σ 2L
σ 2S

73.8 ± 9.4
7.38 ± 0.070−0.176 ± 0.031
4.80
0.238
0.215

46.4 ± 6.4
7.23 ± 0.091
0.169 ± 0.042
1.1
0.281
0.253

24.3 ± 1.10
6.80 ± 0.031−0.107 ± 0.091
0.25
0.289
0.277

SF > 1.1 mm/yr
SF ≤ 1.1 mm/yr

6
3L
5
1
10

MC (

= 0.28110=10.6 bars, CLD =1.4, Wmax =30) [SF =1.1 mm/yr]
2

10

3

LE , km1
Residual from Top Plot
M w = 0.144 log(SF / 1.1 mm/yr )

Mw

0.5
0
–0.5
–1 3S
10 –3

= 0.255
10 –2

10 –1

10 0

10 1

10 2

Model M1: The Linear Model

S , mm/yr
F

(c)

9

NormalMw

8

7SF > 0.25 mm/yr
SF ≤ 0.25 mm/yr

6
3L
5
1
10

MC (

= 0.31210=14 bars, CLD =1.2, W max =18) [SF =0.25 mm/yr]
2

10

3

10

1

LE , kmResidual from Top Plot
M w = –0.056 log(SF / 0.25 mm/yr )

1

M

w

0.5
0
–0.5
–1 3S
–2
10

= 0.305
10

–1

10

0

S , mm/yr
F

Figure 7.

Model M3 (equation 3) for (a) strike-slip, (b) reverse,and (c) normal faults. Coefficients for the lines are given in Table 5.
Other figure details are the same as in Figure 5.

lower frame shows the residuals from this prediction,
defined as
EQ-TARGET;temp:intralink-;df7;55;149^ wi ;
δM wi  M wi − M7

for each considered earthquake, and the solid line is given byδMw  c2 logSF =S0 . Model coefficients and uncertainties
in estimates of Mw for models M1–M3 are given in
Tables 3–5, respectively.

The parameters for the linear models are given in
Table 3. Figure 8 shows the distribution of coefficients foundfor 10,000 trials for strike-slip faults. The widths of these
distributions are used to estimate the uncertainty in eachcoefficient. The coefficients c0 , c1 , and c2 are found simultaneously, as opposed to a possible alternative approach, inwhich c0 and c1 could be found first, and then c2 is determined by a second independent linear fit to the residuals.
For strike-slip events, which dominate the data,c2  −0:198  0:023 (Fig. 5a) so δM w is observed to be adecreasing function of slip rate, similar to AWS96. The datawith a reverse mechanism support δM w increasing, rather
than decreasing, with increased slip rate (Fig. 5b), whereas
for the events with a normal mechanism the slip-rate dependence of δM w is not distinguishable from zero (Fig. 5c).
Considering the distribution of slip-rate data for reverse
faults in Figure 5b, it may be observed that the finding ofslip-rate dependence is the result of mainly a single outlier,the Marryat earthquake (M w 5.8, event number 28 in Table 1)
which is reported to have a slip rate of 0:005 mm=yr. Intracontinental events are included considering, based on
Byerlys law, that the physics of rupture of crystalline rocks
within the range of typical crustal compositions is not,a priori, different merely because the fault is located far froma plate boundary or that rock type might be different (e.g.,
Byerlee, 1978; Scholz, 2002). Also, the Marryat Creek event
tends to decrease the slip-rate dependence of δM w , as aconsideration of the remaining points would reveal.

−119−

2570
500

J. G. Anderson, G. P. Biasi, and S. G. Wesnousky

Strike slip

0
4.5

4.55

4.6

4.65

4.7

4.75

4.8

4.85

4.9

4.95

5

c0
500
0
1.15

1.2

1.25

1.3

1.35

1.4

1.45

c1

strike-slip case plus some noise will still resolve to a decently
significant trend even though we added only noise. In applying the combined regression to dip-slip faults, we may be projecting back from the strong case into the noise, and saying
things about future dip-slip earthquake expectations that are
not likely based on the available data.

500
0
–0.3

–0.28

–0.26

–0.24

–0.22

–0.2

–0.18

–0.16

–0.14

–0.12

500
0
0.15

0.2

0.25

0.3

0.2

0.25

0.3

σ 1L

0.35

0.4

0.45

0.5

0.35

0.4

0.45

0.5

500
0
0.15

σ 1S

Model M2: The Bilinear Model

–0.1

c2

Figure 8. Coefficient distributions for the linear strike-slip
model (M1, equation 1). The bar chart shows number of occurrences of parameter values among 10,000 realizations for randomly
selected values of Mw , LE , and SF within the uncertainty range of
each. The solid gray line shows the mean value of each parameter.
The dashed gray line shows the value found for the preferred valueof Mw , LE , and SF for each earthquake. The clear negative value ofc2 corresponds to decreasing relative magnitude predictions withincreasing slip rate.

Table 4 gives estimated coefficients for model M2(equation 2), and Figure 6 illustrates the fit to the data. Without the slip-rate adjustment, the bilinear model fits theobserved magnitudes as well or better than the linear model
M1, as shown by similar or smaller values of σ 2L than thecorresponding values of σ 1L . The results again show adependence of magnitude on slip rate for strike-slip faults(Fig. 6a) but not dip-slip faults (Fig. 6b,c). The value of σ 2S
is comparable to σ 1S for the strike-slip case but larger for the
dip-slip faults. For the strike-slip case, the fit to the data in
Figure 6a is better at large rupture lengths than in Figure 5a.
Model M3: The Constant Stress-Drop Model

Parameters for model M3 are given in Table 5, and the fit
to the data is illustrated in Figure 7. Some features of Figure 7
are noteworthy. For the strike-slip case, the points for faults
with low slip rates (solid points) are mostly above the averNonetheless, the positive slope of δM w in Figure 5b forage model, whereas points with high slip rates (open circles)
reverse faulting is not a robust result.
are mostly below the average. This slip-rate dependence isConsidering Figure 5, the results for the linear model proreinforced in the lower frame of Figure 7a, where the slope
vide an indication that it is not appropriate to combine differof the linear fit to the residuals is more than five standard
ent fault mechanisms in this type of regressions. The AWS96
deviations of the slope different from zero. The variance remodel from all rupture mechanisms was Mw  5:12
duction by the addition of the slip-rate term is statistically
1:16 log LE − 0:20 log SF , which is only slightly different
significant with 80% confidence, based on the F-test. The
from the strike-slip case in Figure 5a. That result is consistentsame remarks apply for models M1 (Fig. 5) and M2 (Fig. 6).
with the AWS96 model being dominated by strike-slip earthThe strike-slip case in Figure 7a uses W max  15 km,
quakes, and thus demonstrates continuity with the previous
whereas Table 5 gives model M3 parameters forstudy. However, results here separated by mechanism indicate
W max  20 km as well. The data do not prefer either of thesethat the slip-rate dependence in AWS96 is controlled by thetwo models, as the curves and the misfits characterized bybehavior of strike-slip earthquakes, and not much affected by
σ 3L and σ 3S are barely distinguishable, so the plot for thethe normal mechanisms that show little or no slip-rate dependW max  20 km model is not shown. For the 20 km wide
ence, and the reverse mechanisms that potentially show astrike-slip case, both σ 3L and σ 3S are smaller than the equivdifferent dependence. Suppose as a thought experiment thatalent uncertainties in models M1 or M2. Although this imthe dip-slip cases have no slip-rate dependence, or in otherprovement is small and statistically insignificant, it iswords, that the variability with slip rate is pure noise. A strong
encouraging that a model with constantstress drop achieves this result. Hanks
Table 5and Bakun (2014) discussed the difficulCoefficients for Model M3, Use in Equation (3), the Different Fault Types
ties associated with several scaling models
Considered Separately, and Earthquakes Listed in Tables 1 and 2
for long strike-slip faults, which fit the lonParameter
Strike Slip (15 km)
Strike Slip (20 km)
ReverseNormalgest earthquakes either by increasing the
rupture width by penetrating into the crustΔτC (bars)
24.9 ± 1.1
15.3 ± 0.7
10.6 ± 0.7
14.0 ± 1.5
below the depths of microearthquakes or
CLW
3.8
2.9
1.4
1.2
W max (km)
15
20
30
18
by increasing the stress drop. Althoughc2−0.170 ± 0.029 −0.174 ± 0.029 0.144 ± 0.027 −0.056 ± 0.095
Hanks and Bakun (2014) consider theS0 (mm=yr)
4.8
4.8
1.1
0.25
deep penetration of strike-slip faults below
σ 3L
0.236
0.235
0.281
0.312
the depth of microearthquakes to be
σ 3S
0.214
0.210
0.255
0.305
unlikely, we provide both models. Recent−120−Fault-Scaling Relationships Depend on the Average Fault-Slip Rate(a) 8.5
8

Strike slip, S F = S 0 = 4.8 mm/yr

Mw

7.5
7
M1
M2
M3

6.5
6
5.5
101

(b) 8.5
8

10 2

10 3

Reverse, S F = S 0 = 1.1 mm/yr

Mw

7.5

this corresponds to ΔτS ≈ 50 bars (see the Appendix) in the
more frequently used model of Kanamori and Anderson(1975). We suggest that variability of the aspect ratio must
contribute to the uncertainties in these scaling relations at thelower magnitudes.
For the reverse-faulting data, we considered values ofW max up to 30 km because reverse faults can have low dips,and that upper limit is the preferred value. For normal faulting,
we only considered the values of W max > 18 km, because
constrained observations of normal faults imply that the fault
width can be that wide (e.g., Richins et al., 1987). For modelM3 to fit the sparse normal-faulting data as well as model M2,
we would need to use a much smaller value of W max .

7
6
5.5
101

(c) 8.5
8Comparisons

M1
M2
M3

6.5

10 2

10 3

Figure 9 compares the models for the three different
types of mechanisms. Models M2 and M3 tend to resemble
each other most closely, whereas model M1, being linear,tends to give larger magnitudes for long and short faults
but smaller magnitudes in the center of the length range.

Normal, S = S = 0.25 mm/yr
F

0

7.5

Mw

2571

Discussion

7
M1
M2
M3

6.5
6
5.5
1
10

10

2

10

3

Rupture Length, km
Figure 9. Comparisons of models M1, M2, and M3 for(a) strike-slip, (b) reverse, and (c) normal faults, using SF  S0
for each mechanism.
studies that favor deep penetration include Graves and Pitarka (2015) based on experience in modeling ground motions near the fault and Jiang and Lapusta (2016) basedon the seismic quiescence of the ruptures of past large earthquake such as the 1857 earthquake on the San Andreas fault.
The ability to model the data using equation (3) unfortunately does not resolve the “no high stress drop/no deep slip
enigma” articulated by Hanks and Bakun (2014), but rather
pushes it into issues with the aspect ratio and the absolutevalue of the constant stress drop. The W max  15 km model
uses a large aspect ratio of CLW  3:8, compared with
CLW  2:9 for the W max  20 km model, or 2.4 at thetransition to fix the width of 15 km in the Hanks and Bakun(2014) model. The higher aspect ratio for the W max  15 km
model would also imply that earthquakes such as the M w 6.6
Superstition Hills event (number 25) or the M w 6.2 Parkfield
1966 event (number 49) only penetrate from the surface
to about 7 km depth. Also, the stress drop for the W max 
15 km model is rather high, ΔτC ≈ 25 bars, considering that

The larger data set modeled here compared with AWS96
expands our understanding of slip-rate dependence for thescaling of magnitude and rupture length. Improvements inestimates of the magnitudes of earthquakes are realized withslip-rate dependence for strike-slip faults for all three models
considered here. Thus, the slip-rate dependence in this case is
not an artifact of the underlying scaling model. The distribution of data in LE − SF space (Fig. 4b) gives further reason
for confidence in the strike-slip case. On the other hand, individual models for reverse and normal faulting have, at best,
an equivocal place for slip-rate dependence. This again isconsistent with the uneven distribution of data on the plotsof LE − SF in Figure 4c,d. For models M1 and M2 of thenormal-faulting events, but not for model M3, the sign ofslip-rate dependence agrees with the strike-slip case. Thus,
normal faulting could have slip-rate dependence nudging estimates toward smaller magnitudes for higher slip-rate faults
but lacks sufficient data to prove it. The reverse-faulting
events disagree even at the sign of the effect. The disagreement is present whether we retain either or both of the apparent outliers in Figures 5–7. Thus, based on the current
data, we do not find support for the general reduction of magnitude with slip rate implied by the combined set regressionof ASW96. It appears that the strength of the slip-rate effectamong strike-slip events and their sheer numbers relative to
dip-slip events overwhelm the ambiguous (normal) and contrary (reverse) data, leading to an apparently general slip-rate
relationship among all data. Thus, our new data set contributes to the understanding that slip-rate dependence is dominantly a strike-slip fault effect that is not inconsistent withnormal faulting, and not apparently consistent with reverse

−121−

2572

J. G. Anderson, G. P. Biasi, and S. G. Wesnousky

mechanism fault rupture. The data available to AWS96 did
not permit this distinction.
If we are guided by studies of earthquake source
physics, model M3 may be preferable to models M1 or M2.
Specifically, the advantage would be the constant stress dropof earthquakes over the full range of magnitudes, consistent
with, for example, Allmann and Shearer (2009). The slope ofthe linear model M1 with rupture length implies that stressdrop increases significantly with rupture length for large
earthquakes. The slopes of the bilinear model M2 are consistent with simple models for scaling with constant stressdrop in the small and large earthquake domains, but the stressdrops in the two domains are different. In addition, because
the buried circular rupture model by construction does notreach the surface, its applicability to the short ruptures ofmodel M2 is not obvious.
Constant stress-drop model M3 has the important advantage compared with dislocation models in an unbounded
space in that it is explicitly designed for surface-rupturing
earthquakes. For this reason, we might expect that it will perform well where magnitude scaling is required for calculation of synthetic ground motions (e.g., Goulet et al., 2015).
The Chinnery (1964) model has uniform slip with a singularity of stress drop near its edges, which enables a closedform solution. Stress drop in actual earthquakes is a variable
function of location on the fault, so single values are always
averages. The Chinnery (1964) approach is probably as reasonable as others. The application of the same functional
form for dip-slip earthquakes is entirely ad hoc, of course.
Although it is more complicated, its consistency with aphysical model with a constant stress drop commends it
as a preferred regression. Compared with the better-known
equations summarized by Kanamori and Anderson (1975),
the stress-drop parameter in this model is smaller, emphasizing that average stress-drop estimates are model dependent.
The adjustment that decreases magnitude for high-sliprate strike-slip faults implies that the stress drop on thosefaults is lower than on faults of the same length with lower
slip rate. The finding is consistent with the observations ofKanamori and Allen (1986) and Scholz et al. (1986) that alonger healing time results in a larger stress required to initiate rupture and thus a higher stress drop. For normal or reverse faulting, the slip-rate dependence is low, and the sliprate coefficient c2 is indistinguishable from zero. The findings suggest that, if c2 is not zero for these cases, then c2 is
positive for reverse-faulting earthquakes. This is contrary to
the hypothesis of Kanamori and Allen (1986). We suggestthat if this positive slope is confirmed with added data, thephysical mechanism may be related to the dynamics of rupture. For a reverse fault, the dynamic stresses on a rupture
propagating up-dip are tensile as rupture approaches the surface, so the coefficient of friction or cohesion on the fault is
less relevant.
There are a number of future studies that should be performed to improve upon the results presented here. The firstis to examine the consistency of the models, and especiallymodel M3, with the observed fault displacement. If the results, based on the definition of seismic moment (equation A2), agree with seismic data, the scaling relationshippresented here would be an alternative to the self-consistentscaling model of Leonard (2010, 2014) for earthquakes incontinental crust.
The second issue that deserves attention is handling
multisegment faults. We consider, for instance, the 1905 Bulnay, Mongolia, earthquake, which is the strike-slip point in
Figures 5–7 at 375 km and M w 8.5. The 375 km length is the
distance from one end of the rupture to the other, and does
not include a spur fault in between that is 100 km long. This
event points out that the standard deviations σ xL and σ xS for
all three models include the potential presence of spur or subparallel faults that do not increase the total end-to-end lengthof the rupture. Several other faults in Tables 1 and 2 have
similar issues. A better understanding of how seismic moment is distributed on multiple segments and fault splays,
as well as how best to measure the lengths of multiple segment ruptures and how to recognize these features ahead ofthe earthquake would help to reduce uncertainties in futurestudies of scaling relations. If the result is different from theapproach used by Uniform California earthquake rupture
forecast, v. 3 (UCERF3), it could have a direct impact on
future seismic-hazard analyses.

Conclusions
The primary question asked by this research is if the introduction of slip rate on a fault helps to reduce the uncertainties in estimates of magnitude from observations ofrupture length. We find that such a slip-rate dependence isreasonably well established for strike-slip cases: as the slip
rate increases for any given fault length, the predicted magnitude tends to decrease. This result is robust in the sense thatthe slope of the residuals with slip rate is significantly different from zero, and the variance reduction is modestly significant for all three of the considered models relating rupture
length and magnitude. For reverse- and normal-faulting
mechanisms, on the other hand, our data do not demonstrate
the presence of a significant slip-rate effect in the relationshipbetween rupture length and magnitude. Compared withoriginal results in AWS96, we now suggest slip rate be included only for strike-slip faults.
The constant stress-drop model presented here has potential for progress on a standing difficulty in ground-motion
modeling of an internally consistent scaling of magnitude,
length, down-dip width, and fault displacement. Current relations in which magnitude scales with length or area lead tounphysical stress drops or unobserved down-dip widths, respectively. By working from the model of Chinnery (1964),
our constant stress-drop model has the advantage of startingwith realistic physics including the stress effects of surface
rupture. Work remains to be done in comparing displacements predicted from our model with observations, but
the fact that it fits the current magnitude-length-slip-rate data

−122−

Fault-Scaling Relationships Depend on the Average Fault-Slip Rateas well or a bit better than the linear and bilinear models suggests that the constant stress-drop model is preferable tomodels that do not have this feature.

Data and Resources
The article “Total least squares” is available at https://en.
wikipedia.org/wiki/Total_least_squares (last accessed February 2015).

Acknowledgments
The authors thank James Brune for many helpful discussions on thisresearch project. The authors thank Tom Hanks, Ivan Wong, and two anonymous reviewers for very helpful reviews. This research was supported in partby the Southern California Earthquake Center (SCEC, Contribution Number
7485). SCEC is funded by National Science Foundation (NSF) Cooperative
Agreement EAR-1033462 and U.S. Geological Survey (USGS) Cooperative
Agreement G12AC20038. This research was supported primarily by the
USGS, Department of the Interior, under USGS Award Number
G14AP00030.
The views and conclusions contained in this document are those of theauthors and should not be interpreted as necessarily representing the official
policies, either expressed or implied, of the U.S. Government.

References
Abercrombie, R. (1995). Earthquake source scaling relationships from -1 to
5 ML using seismograms recorded at 2.5-km depth, J. Geophys. Res.
100, 24,015–24,036.
Allen, C. R. (1975). Geological criteria for evaluating seismicity, Bull. Geol.
Soc. Am. 86, 1041–1057.
Allmann, B. P., and P. M. Shearer (2009). Global variations of stress drop formoderate to large earthquakes, J. Geophys. Res. 114, no. B01310, doi:
10.1029/2008JB005821.
Anderson, J. G. (1986). Implication of attenuation for studies of the earthquake source, Earthquake Source Mechanics, Maurice Ewing Series 6,
Geophysical Monograph 37, American Geophysical Union, Washington, D.C., 311–318.
Anderson, J. G., S. G. Wesnousky, and M. Stirling (1996). Earthquake size
as a function of fault slip rate, Bull. Seismol. Soc. Am. 86, 683–690.
Baltay, A., S. Ide, G. Prieto, and G. Beroza (2011). Variability in earthquake
stress drop and apparent stress, Geophys. Res. Lett. 38, L06303, doi:
10.1029/2011GL046698.
Baltay, A., G. Prieto, and G. C. Beroza (2010). Radiated seismic energy fromcoda measurements and no scaling in apparent stress with seismic moment, J. Geophys. Res. 115, no. B08314, doi: 10.1029/2009JB006736.
Biasi, G. P., and S. G. Wesnousky (2016). Steps and gaps in ground ruptures:
Empirical bounds on rupture propagation, Bull. Seismol. Soc. Am. 106,
1110–1124, doi: 10.1785/0120150175.
Blaser, L., F. Krueger, M. Ohrnberger, and F. Scherbaum (2010). Scaling
relations of earthquake source parameter estimates with special focus
on subduction environment, Bull. Seismol. Soc. Am. 100, 2914–2926.
Bormann, P., and D. Di Giacomo (2011). The moment magnitude Mw andthe energy magnitude Mc : Common roots and differences, J. Seismol.
15, 411–427.
Bormann, P., M. Baumbach, G. Bock, H. Grosser, G. L. Choy, and J. Boatwright (2005). Seismic sources and source parameters, in New Manualof Seismological Observatory Practice (NMSOP), P. Bormann (Editor), Deutsches GeoForschungsZentrum GFZ, Potsdam, Germany,
1–94.
Bormann, P., S. Wendt, and D. Di Giacomo (2013). Seismic sources andsource parameters, in New Manual of Seismological Observatory
Practice 2 (NMSOP2), P. Bormann (Editor), Chapter 3, Deutsches
GeoForschungsZentrum GFZ, Potsdam, Germany, 1–259.

2573

Byerlee, J. D. (1978). Friction of rocks, Pure Appl. Geophys. 116, 615–626.
Castellaro, S., and P. Bormann (2007). Performance of different regression
procedures on the magnitude conversion problem, Bull. Seismol. Soc.
Am. 97, 1167–1175, doi: 10.1785/0120060102.
Castellaro, S., F. Mulargia, and Y. Y. Kagan (2006). Regression problems for
magnitudes, Geophys. J. Int. 165, 913–930, doi: 10.1111/j.1365246X.2006.02955.x.
Chinnery, M. A. (1963). The stress changes that accompany strike-slip faulting, Bull. Seismol. Soc. Am. 53, 921–932.
Chinnery, M. A. (1964). The strength of the Earth’s crust under horizontal
shear stress, J. Geophys. Res. 59, 2085–2089.
Dieterich, J. H. (1972). Time dependent friction in rocks, J. Geophys. Res.
20, 3690–3704.
Goulet, C. A., N. A. Abrahamson, P. G. Somerville, and K. E. Wooddell(2015). The SCEC broadband platform validation exercise: Methodology for code validation in the context of seismic-hazard analyses,
Seismol. Res. Lett. 86, 17–26, doi: 10.1785/0220140104.
Graves, R., and A. Pitarka (2015). Refinements to the Graves and Pitarka(2010) broadband ground-motion simulation method, Seismol. Res.
Lett. 86, 75–80, doi: 10.1785/0220140101.
Hanks, T. C., and W. H. Bakun (2002). A bilinear source-scaling model for
M-logA observations of continental earthquakes, Bull. Seismol. Soc.
Am. 92, 1841–1846.
Hanks, T. C., and W. H. Bakun (2008). M-logA observations for recent large
earthquakes, Bull. Seismol. Soc. Am. 98, 490–494.
Hanks, T. C., and W. H. Bakun (2014). M-logA models and other curiosities,
Bull. Seismol. Soc. Am. 104, 2604–2610, doi: 10.1785/0120130163.
Hanks, T. C., and H. Kanamori (1979). A moment magnitude scale, J. Geophys. Res. 84, 2348–2350.
Hillers, G., and S. G. Wesnousky (2008). Scaling relations of strike-slip
earthquakes with different slip-rate-dependent properties at depth,
Bull. Seismol. Soc. Am. 98, 1085–1101, doi: 10.1785/0120070200.
Ide, S., G. C. Beroza, S. G. Prejean, and W. L. Ellsworth (2003). Apparent
break in earthquake scaling due to path and site effects on deep borehole
recordings, J. Geophys. Res. 108, 2271, doi: 10.1029/2001JB001617.
Iida, K. (1959). Earthquake energy and earthquake fault, Nagoya University,
J. Earth Sci. 7, 98–107.
Jiang, J., and N. Lapusta (2016). Deeper penetration of large earthquakes onseismically quiescent faults, Science 352, 1293–1297.
Kanamori, H. (1977). The energy release in great earthquakes, J. Geophys.
Res. 82, 2981–2987.
Kanamori, H., and C. R. Allen (1986). Earthquake repeat time and average
stress drop, in Earthquake Source Mechanics, S. Das, J. Boatwright,
and C. H. Scholz (Editors), Geophysical Monograph 37, 227–235.
Kanamori, H., and D. L. Anderson (1975). Theoretical basis of some empirical relations in seismology, Bull. Seismol. Soc. Am 65, 1073–1095.
King, G. C. P., and S. G. Wesnousky (2007). Scaling of fault parameters forcontinental strike-slip earthquakes, Bull. Seismol. Soc. Am. 97, 1833–
1840.
Leonard, M. (2010). Earthquake fault scaling: Self-consistent relating of
rupture length, width, average displacement, and moment release, Bull.
Seismol. Soc. Am. 100, 1971–1988, doi: 10.1785/0120090189.
Leonard, M. (2012). Erratum to “Earthquake fault scaling: Self-consistent
relating of rupture length, width, average displacement, and moment
release,” Bull. Seismol. Soc. Am. 102, 2797.
Leonard, M. (2014). Self-consistent earthquake fault-scaling relations: Update and extension to stable continental strike-slip faults, Bull. Seismol.
Soc. Am. 104, 2953–2965, doi: 10.1785/0120140087.
Perrin, C., I. Manighetti, J.-P. Ampuero, F. Cappa, and Y. Gaudemer (2016).
Location of largest earthquake slip and fast rupture controlled byalong-strike change in fault structural maturity due to fault growth,
J. Geophys. Res. 121, 3666–3685, doi: 10.1002/2015JB012671.
Richins, W. D., J. C. Pechmann, R. B. Smith, C. J. Langer, S. K. Goter, J. E.
Zollweg, and J. J. King (1987). The 1983 Borah Peak, Idaho, earthquake and its aftershocks, Bull. Seismol. Soc. Am. 77, 694–723.
Rolandone, F., R. Burgmann, and R. M. Nadeau (2004). The evolution of theseismic-aseismic transition during the earthquake cycle: Constraints

−123−

2574

J. G. Anderson, G. P. Biasi, and S. G. Wesnousky

from the time-dependent depth distribution of aftershocks, Geophys.
Res. Lett. 31, L23610, doi: 10.1029/2004GL021379.
Sato, R. (1972). Stress drop for a finite fault, J. Phys. Earth 20, 397–407.
Scholz, C. H. (1982). Scaling laws for large earthquakes: Consequences forphysical models, Bull. Seismol. Soc. Am. 72, 1–14.
Scholz, C. H. (2002). The Mechanics of Earthquakes and Faulting, Second
Ed., Cambridge University Press, Cambridge, United Kingdom, 471 pp.
Scholz, C. H., C. A. Aviles, and S. G. Wesnousky (1986). Scalingdifferences between large intraplate and interplate earthquakes, Bull.
Seismol. Soc. Am. 76, 65–70.
Shaw, B. E. (2009). Constant stress drop from small to great earthquakes inmagnitude-area scaling, Bull. Seismol. Soc. Am. 99, 871–875.
Shaw, B. E., and C. H. Scholz (2001). Slip-length scaling in large earthquakes: Observations and theory and implications for earthquake physics, Geophys. Res. Lett. 28, 2995–2998.
Shaw, B. E., and S. G. Wesnousky (2008). Slip-length scaling in large earthquakes: The role of deep-penetrating slip below the seismogenic layer,
Bull. Seismol. Soc. Am. 98, 1633–1641.
Strasser, F. O., M. C. Arango, and J. J. Bommer (2010). Scaling of the source
dimensions of interface and intraslab subduction-zone earthquakes
with moment magnitude, Seismol. Res. Lett. 81, no. 6, 941–950.
Tocher, D. (1958). Earthquake energy and ground breakage, Bull. Seismol.
Soc. Am. 48, 147–153.
Wells, D. L., and K. J. Coppersmith (1994). New empirical relationships
among magnitude, rupture length, rupture width, rupture area, and surface displacement, Bull. Seismol. Soc. Am. 84, 974–1002.
Wesnousky, S. G. (1988). Seismological and structural evolution of strikeslip faults, Nature 335, 340–343.
Wesnousky, S. G. (1999). Crustal deformation processes and the stability of the
Gutenberg-Richter relationship, Bull. Seismol. Soc. Am. 89, 1131–1137.
Wesnousky, S. G. (2008). Displacement and geometrical characteristics ofearthquake surface ruptures: Issues and implications for seismic-hazard analysis and the process of earthquake rupture, Bull. Seismol. Soc.
Am. 98, 1609–1632.
Wesnousky, S. G., C. H. Scholz, and K. Shimazaki (1983). Earthquake frequency distribution and the mechanics of faulting, J. Geophys. Res. 88,
9331–9340.

Appendix

EQ-TARGET;temp:intralink-;dfa2;313;733

 E  μLE W E D E;
M 0  μAE D

A2

in which μ is the shear modulus, AE is the fault area ruptured
 E is the average slip over that area.
in the earthquake, and D
For a fault that is approximately rectangular AE  LE W E, inwhich LE is the rupture length measured along strike and W E
is the down-dip rupture width.
Substituting equation (A2) into (A1), one obtains (forcgs units)

EQ-TARGET;temp:intralink-;dfa3;313;626222
 E  2 log μ − 16:1:M w  log LE  log W E  log D
3
3
3
3

A3

This justifies the models that relate magnitude to the log offault length, width, and mean slip. Slopes different from 2/3
result from correlations among the fault parameters LE , W E , E . Wells and Coppersmith (1994) found that the modeland D
EQ-TARGET;temp:intralink-;dfa4;313;520

M w  c1 log LE  c0

A4

predicts magnitude from rupture length with a standard
deviation of the misfit σ 1  0:28.
The possible dependence of stress drop or magnitude on
slip rate was recognized by Kanamori and Allen (1986) and
Scholz et al. (1986). With the addition of the slip-rate term,
equation (A4) becomes, used by Anderson et al. (1996; hereafter, AWS96):
EQ-TARGET;temp:intralink-;dfa5;313;403

Mw  c0  c1 log LE  c2 log SF :

A5

Testing for a logarithmic dependence on the geological faultslip rate SF can be motivated by findings in Dieterich (1972).
In this article, equation (A5) is equivalent to model M1.

Fault-Scaling Relations

Constant Stress-Drop Scaling

BasicsThis article proposes models to estimate the moment
magnitude of earthquakes based on observed surface rupture
lengths and slip rates. The moment magnitude definition that
we use is implicit in Kanamori (1977):2
log M0 − 16:1 :
Mw 
A1
3
EQ-TARGET;temp:intralink-;dfa1;55;231

The units of seismic moment M0 are dyn·cm in equation (A1). This definition differs slightly from the equation
used by Hanks and Kanamori (1979) but is the equation
recommended for seismic network operations by theInternational Association of Seismology and Physics ofthe Earth’s Interior since 2005 (see Bormann et al., 2005,
2013; Bormann and Di Giacomo, 2011; and references
therein), and thus is the relationship used by most seismic
networks throughout the world. The seismic moment is defined as

The static stress drop ΔτS is the average decrease in theshear stress acting on the fault as a result of the earthquake
and is proportional to the ratio of average slip to a fault dimension. Seismic observations have found that the averagevalue of ΔτS is approximately constant (∼4 MPa, ∼40 bars)
over a broad range of earthquake magnitudes (e.g., Kanamori
and Anderson, 1975; Allmann and Shearer, 2009), although
there is considerable scatter in these data. Seismic moment,
and thus Mw through equation (A1), can be expressed as afunction of fault dimension and stress drop, as recognized by
Kanamori and Anderson (1975). Selected models are summarized in Table A1.
The equations in Table A1 indicate that constant stressdrop implies the slope c1  2:0 for small faults (first case 1)
when LE is equated to the diameter of the circular fault andc1  2=3 for long faults (second and third cases). These observations motivate a bilinear approach to fit the data, which
is model M2 in this article. The bilinear approach is formulated as follows:−124−Fault-Scaling Relationships Depend on the Average Fault-Slip Rate

2575

Table A1

motivates the development of the modelthat is described in the next section.

Models from Kanamori and Anderson (1975) for the Relationship of Fault Size,Stress Drop, and M w
CaseBuried, circularM016
3
7 Δτ S RE

Relations Based on Chinnery (1964)

Implied Magnitude Relations*Mw  log AE  23 log ΔτS  3:0089
If LE  2RE : Mw  2 log LE  23 log ΔτS  2:904
Mw  23 log LE  43 log W E  23 log ΔτS  3:1359
Mw  23 log LE  43 log W E  23 log ΔτS  3:3141

Chinnery (1963, 1964) calculated a
stress drop for a rectangular strike-slip
Strike slip, long
fault that ruptures the surface. Unlike
Dip slip, longthe circular slip model, the free surface*AE  πR2E , fault area in km2 ; RE , fault radius; W E , fault width; and LE , fault length in km, andin the Chinnery model is present for smallΔτS , stress drop in bars.
earthquakes. His equations assume a uniform slip on the fault. Thus, the stress drop
 
 
L
S
is variable over the fault and becomes singular at the edge ofLE < Lbp
Mw  M bp c1C log E c2 log F
Lbp
S0
the fault. His equations give the stress drop on the surface at  
the midpoint of the rupture. Numerical solutions in ChinneryLESFc2 log
LE ≥ Lbp ; A6
Mw  M bp c1L log(1963) show relatively uniform stress drop over large porLbp
S0
tions of the fault. Chinnery (1963) thus suggests that the results are valid to represent the fault stress drop so longin which the length Lbp is the length at which the length
because the zone of slip fall-off is much smaller than the area
dependence of the scaling relationship changes from theof the fault. The key advantage provided by this approach issmall fault model with slope c1C  2 to the long fault model
to provide a useful analytical solution.
with slope c1L  2=3. The slip rate S0 is a reference slip rate
For the rectangular fault with length LE , width W E , andwhich can be chosen arbitrarily, but is conveniently chosen to
aspect ratio CLW  LE =W E , the stress drop in the Chinnery
be the log average slip rate in the data, so that setting
model ΔτC at the midpoint at the surface isSF  S0 gives the best fit when slip rate is unknown. The
E
μD
constant M bp is the magnitude corresponding to a fault withC L ; W ;
ΔτC 
A7
length LE  Lbp and slip rate SF  S0 . Equation (A6) has
2π 1 h E
three unknown coefficients (Mbp , Lbp , and c2 ), which isin whichthe same number as in equation (A5).

2Lh
3 Lh 3a  4W E 
However, there are issues with the applicability of the: A8
C1 Lh ; W E  −equations in Table A1. The foremost, for the long faults, isaW E Lhaa  W E 2
the width of the seismogenic zone. Table A1 shows that W E
Note that Lh  LE =2 and a  L2h  W 2E 1=2 . Observe thatis twice as influential as the fault length, so it needs to be
C1 has dimensions of 1=length, and thus C−1
considered carefully. One approach to estimate this width
1 is effectively
thefault dimension that is used for calculating the strain. In
is to use the maximum depth of microearthquakes. By this E C1 .
other words, the strain change in the earthquake is ∼D
approach, for strike-slip earthquakes the maximum depth of
An equation for the seismic moment can be obtained by solvmicroearthquakes equates directly to an estimate of the fault
 E and substituting in equation (A2).
ing equation (A7) for D
width, whereas for a reverse or normal fault the dip is incorThe result is
porated. The problem is that the maximum depth of seismogenic rupture in large earthquakes is difficult to observe.
LE W EM 0  2πΔτCA9King and Wesnousky (2007) discuss this difficulty andC1 Lh ; W E 
present arguments for why the down-dip width might beand thus
larger in large earthquakes, at least up to some limit greater 
than that inferred from the depth range of small earthquakes,
2
2
2
2πW E
2
logL
logΔτ
log
16:1:
−
M



w
E
C
because rocks below the depths of mircoearthquakes might3
3
3 C1 Lh ;W E 
3experience brittle failure under high strain rates. If the widthA10
increases in general for long ruptures, stress drop is no longer
as high for these events because stress drop is inversely proAdditional insight into the geometrical term can be obtained
portional to W E , and furthermore the slope c1 can no longer
by observing that a is the length of the diagonal from the midbe reliably constrained by the models in Table A1. King andpoint of the fault at the surface to either of the bottom corners.
Wesnousky (2007) propose that this explains the proposal by
If the dip of this line is γ, then tan γ  W E =Lh  2=CLW ,Scholz (1982) that slip in large earthquakes is more nearly
Lh  a cos γ, W E  a sin γ, and one can rewrite
proportional to rupture length than to rupture width.
Another issue is that the first equation in Table A1 as1
C1 Lh ; W E  
Cγ;
A11
sumes that the circular fault is confined within the Earth andWEthus neglects free surface effects, while by definition all ofin whichthe events considered in this study rupture the surface. Thisπ
2
2 ΔτS W E LEπλ2μ
24λμ ΔτS W E LE

EQ-TARGET;temp:intralink-;dfa6;55;592

EQ-TARGET;temp:intralink-;dfa7;313;451

EQ-TARGET;temp:intralink-;dfa8;313;408

EQ-TARGET;temp:intralink-;dfa9;313;281

EQ-TARGET;temp:intralink-;dfa10;313;236

EQ-TARGET;temp:intralink-;dfa11;313;120

−125−

2576

Thus, one can rewrite equation (A7) as

EQ-TARGET;temp:intralink-;dfa13;55;686

Cγ Dμ E:
ΔτC 
2π W E

M0 

2π
Δτ L W 2 :Cγ C E E

EQ-TARGET;temp:intralink-;dfa15;55;460

LE< W max :CLW

EQ-TARGET;temp:intralink-;dfa16;55;388

LE≥ W max :CLW

M0  πΔτC LE W 2max :

LD

6

10

1

10

C

=10 bars

C=30 bars

C=100 bars

2

L, km

(b)
8

W

max

=15 km;C=30 bars

7
CLD =1

6

CLD =2
CLD =4

5
4
10

1

10

2

L, km

(c)

8

C =1,
LD

C

=30 bars

7

A15

A17



8
LE
2
2
2π
log

logΔτ

−16:1
2logL
<
E 3
C 3
CLW <W max
C2LW Cγ


Mw 
:2
2πW 2maxLE2
2
3 logLE  3 logΔτC  3 log Cγ −16:1
CLW ≥W max :A18
Equation (A18) will be the third model M3 considered inthis study, with the addition of a slip-rate contribution
c2 logSF =S0 , to the two branches of the equation. The

W max =10

6

W max =15
W max =20

5
4

101

102

L, km

A16

Equation (A17) differs from the second case in Table A1 forthe long strike-slip fault by a factor of 2 (ΔτS  2ΔτC ), in
which the difference is due to the different boundary conditions used for the two solutions at depth.
From equations (A15) and (A16), converting to magnitude, the implied scaling relationship based on the Chinnerymodel isEQ-TARGET;temp:intralink-;dfa18;55;181=15 km; C =1

5

A14

In this case, as the fault length increases while width is held
constant, γ will be decreasing. For the limit of small γ (roughlyγ ≲ 25°), equation (A12) shows that Cγ → 2, so equation (A9) approaches
EQ-TARGET;temp:intralink-;dfa17;55;295

max

7For longer faults, for which the width is limited, equation (A14) becomes
2π
Δτ L W 2
M0 
Cγ C E max

W

4

Because γ, and thus Cγ, depends on the fault aspect ratio,
equations (A9) or (A14), can be used to model a transition
from small-earthquake behavior (e.g., the circular fault inTable A1) to a long-fault behavior. This article, similar toHanks and Bakun (2002), maintains a constant aspect ratio as
the fault length increases, until that aspect ratio implies that thefault width would exceed some maximum. For longer faults,
the width is set to that maximum. Before reaching that maximum, γ and Cγ are constant, and
2π
L3
ΔτC 2E
M0 
Cγ
CLW

8

A13

 E and substituting into equaSolving equation (A13) for Dtion (A2) gives the moment of a vertical strike-slip fault that
ruptures the surface as
EQ-TARGET;temp:intralink-;dfa14;55;607(a)Mw

cos γ sin γ3  4 sin γ: A121  sin γ2

MwCγ  2 cos γ  3 tan γ −Mw

EQ-TARGET;temp:intralink-;dfa12;55;733

J. G. Anderson, G. P. Biasi, and S. G. Wesnousky

Figure A1. Model for Mw based on the Chinnery (1964) scaling as given in equation (A18). (a) Effect of changing the stress dropΔτC . (b) Effect of changing the aspect ratio of the fault. (c) Effect ofchanging the limiting rupture width W max .
unknown parameters in model M3 are ΔτC , CLW , W max ,and c2 . Thus, this model has four parameters to be determined, compared with three parameters in models M1 andM2. Figure A1 shows the effect of the three parametersΔτC , CLW , and W max on magnitude predictions. The stressdrop scales the entire curve upward. The aspect ratio CLD
adjusts the level of the magnitude for short rupture lengths.
The maximum width affects the curvature and how rapidly
the curve approaches the asymptotic slope of 2=3 log LE
for long rupture lengths.
Other Models and Considerations
Sato (1972) overcomes the singularity introduced by
Chinnery (1963, 1964) by assuming a smooth ad hoc slip
function on a finite rectangular/elliptical-shaped fault, andfor that function, calculating the average stress drop resulting
from that slip function. Although the results are informative
for source physics studies, the major disadvantages of this−126−

Fault-Scaling Relationships Depend on the Average Fault-Slip Rate
approach for our application are that the fault is embedded ina whole space, and there is no analytical solution comparable
to equation (A7). Rather, the geometrical factor equivalent toCγ can be computed numerically using equations in Sato
(1972) or read from a figure in the paper. Considering theselimitations, this model was not considered further.
Shaw and Scholz (2001) and Shaw and Wesnousky(2008) implement a numerical model for fault slip in a
half-space with depth-dependent friction. They examine
the statistics of events that rupture the surface. These papers
are interesting for the finding that large surface-rupturing
events also slip below the brittle crustal depths. The scaling
found in the model has properties similar to the scaling in the
Chinnery model. However, the scaling relationship that they
determine has an ad hoc shape, and thus we preferred theanalytical functional form of equation (A18). The physicsbased solution of Chinnery was also preferred to a related
constant stress-drop model by Shaw (2009). This model proposes three regimes of magnitude scaling from length basedon intermediate length–width–displacement–scaling relations and heuristic arguments for transitions between them.
Rolandone et al. (2004) found some empirical evidence
that might be interpreted to support the penetration of rupture

2577

below the brittle seismogenic layer in large earthquakes.
They found that the maximum depth of aftershocks of theLanders earthquake were deeper immediately after themainshock, and then the maximum depth returned to preearthquake levels over the next few years. This might be
explained by high strain rates in the uppermost part of theductile crust, as high strain rates favor brittle failure. However, postseismic strain rates in that depth range would behigh even if seismic rupture of the mainshock did not extend
that deep, so these observations allow, but do not require,dynamic rupture below the long-term average depth of
microearthquakes.

Nevada Seismological Laboratory
University of Nevada
Reno, Nevada 89557
jga@unr.edu
stevew@seismo.unr.edu−127−

Manuscript received 21 November 2016;
Published Online 7 November 2017

Measuring the Improvement in Synthetic
Seismograms Using Shallow Velocity ProﬁlesObtained from Microtremor H/V Spectral
RatiosJohn Anderson and Hiroshi Kawase
October 28, 2018
Abstract
Our ultimate goal is to develop synthetic seismograms that outperformground motion prediction equations through incorporation of the physics
of wave propagation. The present model uses synthetic Green’s functions
generated for ﬂat-layered geological structure. We test four generations ofvelocity models to generate synthetics for the M5.8 Fukushima-Hamadori
earthquake of March 19, 2011. This is a shallow, normal-faulting earthquake that was a foreshock of the M6.7 Fukushima-Hamadori earthquakeof April 11, 2011. This earthquake is strong enough to generate strong
shaking with a high signal-to-noise ratio over a broad frequency band,but small enough to reduce the source complexity. Velocity model VM1 isgiven by JIVSM for the hypocenter of the earthquake. VM2 uses the same
velocity model for all stations, but with improvements to the Q model
for improved spectral matching. VM3 uses instead the JIVSM velocity
model that is given for each station, with the improved Q model. VM4
substitutes a shallow velocity structure based on inversion of H/V frommicrotremors using the inversion method of Kawase et al. (2018). Overall, the quality of the synthetic seismograms, measured as the average
standard deviation of Fourier and pseudoacceleration response spectra,
improves modestly as the models advance from VM1 to VM4. The improvement is marked at some stations, but not at others. A challenge for
future research is to better understand the physical factors behind thesetrends.

1

Introduction

Seismic risk analysis requires hazard estimates in the form of hazard curves
and, increasingly, appropriate seismograms that are compatible with speciﬁc
earthquakes on identiﬁed faults that might aﬀect the site. This paper focuses
on development of realistic synthetic seismograms that can serve as alternatives
to past records.
1

−129−

In the process of developing a model for synthetic seismograms for engineering applications, the ﬁrst question is the selection of an appropriate velocity
model. The SCEC Broadband Platform, as described by Goulet et al. (2015),
used regional models - one for northern California, one for the Mojave Desert,
one for elsewhere in southern California, and two models to represent diﬀerent
regions in Japan. This project evaluates the beneﬁts of going beyond regionalmodels for the purposes of generating synthetic seismograms in two-dimensional
structures.
We went through a number of iterations and model improvements. At theend, we selected four representative models to illustrate the eﬀects of diﬀerencesin the model concept. The result informs the trade-oﬀ between gathering more
information and the reduction of sigma that can result from obtaining that
information.

2

Procedure

2.1

The Earthquake

This case study considers the MW 5.8 earthquake of March 19, 2011 (36.7837
N, 140.5715 E), in the region near the border between Fukushima and Ibaraki
Prefectures . The focal mechanism of this event is normal faulting, on a fault
with a southeast strike (141o ) dipping to the southwest (48o ). Figure 1 shows thelocation of the earthquake and it’s aftershocks in eastern Japan, north of Tokyo.
It was recorded by the Japanese K-NET and KiK-net strong motion network.
Figure 1 locates the 42 stations within 200 km of the fault with estimated valuesof VS30 > 500 m/s, and used in this study.

2.2

Synthetic Seismograms

This project aimed to reproduce key charancteristics of the records at all 42
stations through the use of the composite source model system of generating
synthetic seismograms. This code has been described by Anderson (2015). Synthetics at all frequencies are generated using the representation theorem. Thisrequires a slip model for the source, and Green’s functions to transmit the eﬀectof the source slip to the station. This study focuses on the velocity model used
to calculate the Green’s functions. However, an elementary review of the source
model will be helpful.
The source model consists of the superposition of subevents, where eachsubevent is modeled as the source of a “Brune pulse” (Brune, 1970). Thus eachsubevent has a virtual radius, ri , a time constant τi = ri / (2.34β), and a timefunction Ṁ0 (t) ∼ Δτs te−t/τ . The pulse size is normalized by a subevent stress
drop, Δτs such that the time function integrates to the moment of a circular3
crack of that stress drop: M0i = 16
7 ΔτS ri (Kanamori and Anderson, 1975).
Each subevent has seismic magnitude MW i = 23 (log M0i − 16.1) for momentin dyne-cm (Kanamori, 1977). The numbers and radii of subevents are chosen

2

−130−

randomly from probability distributions that 1) satisfy a statistical GutenbergRichter relationship with b = 1, and 2) sum to the target seismic moment of theearthquake. The largest allowed subevent has rmax = WE /4 where WE is thenarrowest width of the fault. Subevents are placed on the fault at random, but
a circle with radius ri cannot go outside of the fault limits. Each subevent starts
at the time given by the distance from the hypocenter, divided by the rupture
velocity. This model may break down for sites in the near ﬁeld of the large
subevents, but it generally seems to perform well. The total slip represented
by the subevents, found by summing the slip of all the circular subevents inthis study results in variable slip on the fault. From this slip, it is possible toﬁnd the “strong motion generating area”, and verify that it is consistent withthe scaling of the Irikura recipe (Somerville et al., 1999). Furthermore, thatslip distribution is consistent with the slip distribution reported by Kiram et al.
(2016).
The Green’s functions are all generated using a Fortran code ﬁrst written by
Y. Zeng, but subsequently modiﬁed, following the method described by Lucoand Apsel (1983). This method represents the Earth with ﬂat layers, where eachlayer is characterized by the thickness, P-wave and S-wave speeds, attenuation
quality factors Qp and Qs, and density. This study considers the eﬀectivenessof four approaches to determine the velocity model, as summarized in Table 1.
Figures 5 and 6 illustrate the diﬀerences the Q models in VM 1 and VM 2. The
reasons for this change in the Q model will be pointed out in the discussion of
the results.
Table 1: Characteristics of the velocity models considered in this study.
Model
Origin and Main Characteristics
VM 1
Obtained the velocity model for grid point nearest to theepicenter of the chosen earthquake from JIVSM, following
instructions from H. Miyake:
http://www.jishin.go.jp/main/chousa/12_choshuki/dat/nankai/lp2012nankaie_
str.zip.
This model for the crust was used exactly.
VM 2
Velocities are the same as VM1. Q model is modiﬁed from VM1.
Shallow Qp=Qs, as given by the procedure in Anderson (2015).
In the crustal waveguide (depth > 5 km), Qp=Qs=1000.
VM 3
Separate velocity model used for each site. Velocities are fromJIVSM at each station location. Q is determined as in VM2.
VM 4
Like VM3, but shallow velocity structure was replaced withstructure obtained from H/V ratios following the method of
Ducellier et al. (2013) and Kawase et al. (2018).

3

−131−

3

Results

It is useful to begin the discussion of the synthetics with images of selectedsynthetics, compared to data. For this purpose, this paper will focus on thestation IBRH16. As seen on Figure 2, this station is located southwest of thefault, on the hanging wall side of the fault, but beyond the vertical projectionof the fault to the surface. Figures 7, 8, and 9 compare the observed and oneof the model set of synthetic acceleration, velocity, and displacements. For
each station, synthetics were calculated for 50 realizations of the source. The
synthetics for these ﬁgures were computed using VM4. The data and synthetics
have similar amplitudes and durations for all three components of the groundmotion. They also appear to have similar frequency content, although that isdiﬃcult to judge from looking at the time series.
Figure 10 compares the observed Fourier amplitude at station IBRH16 to
the average, over 50 realizations, of the source for VM1, VM2, VM3, and VM4.
Figure 11 shows the equivalent for the pseudo-spectral acceleration (PSA). As aquick visual impression, VM4 appears to have generated the most similar Fourier
spectra, particularly at frequencies near the frequency where the data spectrum
peaks. We attribute the approximate match in frequencies of spectral peaks at
about 4 Hz, 6 Hz, and 8 Hz to the improvements in the velocity model based on
the H/V. VM1 also shows a peak at 6 Hz, which we attribute to coincidence,as VM2 and VM3 are much smaller in this range, and their spectra do notmatch. The VM1 spectrum is too high from ~12-25 Hz, which we attribute to
the Q model. This is one illustration of why the Q model needs to have very
low values near the surface, as in Figure 6 and is used for VM2, VM3, andVM4. The broad trend of the model spectra below 4 Hz is the same, as should
be expected, but none of these models have a convincing match of the ﬁner
structure in this frequency range.
Comparing the observed PSA spectrum with the four models (Figure 11),
the result of VM4 seems to provide the best match. At long periods, all four
spectra converge as expected. At the short period side of the spectrum, VM4
comes closest to predicting the observed peak acceleration. The width andamplitude of the VM$ model peak is closer to matching the data than any of
the others. VM2 and VM3 lack a peak in this period range. VM1 does have apeak that is similar to the data at Tn~0.2 s, but this model also has large peak
between 0.05s and 0.1 s that is not present in the data.
Figure 12 shows the ratio of Data/Model for spectra such as those in Figures
10 and 11. This ﬁgure plots, on a logarithmic axis, the average spectral ratios foreach of the 42 stations in Figure 1 based on the VM4 model as the denominator.
For subsequent analysis, the residual is deﬁned as the ln(data/model). At low
frequencies / long periods, the residual converges to zero because that part ofthe spectrum is controled by the seismic moment of the source. The subevent
stress drop for each velocity model was adjusted to give zero residual for Tn=0.
This is not necessarily the optimal choice, but is useful to for the comparative
purposes of this paper.
One test of the diﬀerent velocity models is to see if there are trends in the
4

−132−

residuals as a function of distance. Perusal of Figure 12 may show some trends,
but it is easier to to see if spectral residuals at a speciﬁc frequency or period
are plotted as a function of distance. This has been done, and Figure 13 shows
the slopes, with their uncertainties, for each model at four frequencies. At 1 Hz
and at ~3 Hz, the slopes are not zero, but a zero slope is within the one-sigma
uncertainty range for all four models. At 10 Hz / 0.1 s, the residuals in the slope
is signiﬁcantly non-zero for VM1, but not for the other three models. Again,
this unacceptable feature has been corrected by the change in the attenuation
structure illustrated in Figures 5 and 6. At the longest periods shown (0.3 Hz/ 3 s), all of the models have a signiﬁcant tendency to have a larger residual atshort distances, and a smaller residual at long distances.
Another way to compare the four velocity models is to compare the mean
spectral residuals. This is shown in Figure 14. In this ﬁgure, the average
response spectral period tends to be better for VM4 than any of the others. In
the residuals for the Fourier spectrum, VM4 shows the largest deviation from
0.3-3 Hz, but tends to be better than the others at other frequencies. The
corresponding standard deviations of the models, derived from the variability ofstation means as in Figure 12, is shown in Figure 15. Here, the diﬀerences seemsmall compared to the overall values. VM4 tends to be best from 1-10 Hz in the
Fourier spectrum, and a corresponding period range in the response spectrum.
Figure 16 attempts to concisely summarize the observations in Figure 15.
The misﬁt is the sum of the absolute values of the residuals in Figure 14, sampled on the log frequency / period axis scale. The standard deviation is theaverage of the standard deviations in Figure 15, sampled at equispaced points
on the log frequency / period axis scale.. The relative values of the standard
deviations indicate that the Fourier amplitude models based on the apropriate
deep structure beneath the site is better than a single regional model, and thatincorporating shallow site information is even better. The impact of including
the shallow site information on the response spectrum does not show subsequent
reduction of sigma.

4

Discussion / Conclusions

There are situations where broadband synthetic seismograms that satisfy thewave equation are useful and ﬂat-layered velocity models are available but three
dimensional models are not yet available. There are also situations where broadband ground motions in ﬂat-layered models are suﬃcient for addressing speciﬁcscientiﬁc questions. This case study is relevant to those situations.
The broadest conclusion is that on average, it seems better to use a velocity
model that is appropriate for the site, including near-surface velocities, if thatis available. In the best case, a set of site-speciﬁc models like those used for
VM4 might, hopefully, provide ground motion estimates with uncertainties thatare comparable to determine single-station sigma in ground motion prediction
equations. Figure 12 shows, however, that the standard deviations found using
the synthetics have standard deviations of their estimated parameters that are

5

−133−

considerably greater than the single-station sigma values found by RodriguezMarek et al. (2011) at KiK-net stations. Consistent with this observation,
review of spectra on a site-by-site basis ﬁnds that the extent to which the models
were improved in the shift from VM3 to VM4 is variable. A future challenge
is to understand why these results are mixed, and consider how they can be
improved to close the gap between the best models of synthetic seismogramsand the results for single-station sigma.
The selected station, IBRH16, was chosen primarily because of it’s location
close to the fault oﬀ the hanging wall. Many urban areas near normal faults
share that geometry. However, as a qualitative impression, in about a third of
the cases, like IBRH16, VM4 was obviously best, in another third the spectral
shapes were changed towards the observations and away from the other models,
but the amplitude of the ﬁt was oﬀ, and in about a third, model VM4 did notstand out from the other models. Indeed, it is clear in Figure 12 that there
are some stations with conspicuously poor agreement with the model. We haveobserved, for instance, that some of the largest residuals are seen at stations
in deep valleys near the west side of Honshu, within the volcanic range. It is
known that Q is relatively low in that region, and perhaps also topography has
an eﬀect. This is one example where the physics aﬀecting the ground motions
is incompletely modeled in this study.
Another shortcoming in the physics included in the ﬂat-layered model isseen in Figure 13, which indicates that the low frequency surface waves are
not attenuating fast enough, on average. Indeed, displacement seismograms
at stations beyond 150 km show strong surface waves that are not present in
the data, and their Fourier spectra show peaks at ~0.2 Hz that are stronger
than the data, It seems likely that scattering by near-surface complexity, thatis not represented in the model, would reduce the amplitude of these surface
waves. The very low Qp and Qs in the near-suface layers apparently does notcause these waves to attenuate strongly enough. This might alternatively becorrected with a frequency-dependent Qp and Qs, but that is not currently
implemented in the Green’s function calculations. There are uncomplicatedpatches to the present model. We have experimented with separate Green’s
functions generated with low Q, and combining them with the high-frequenciesof the models shown here using match ﬁlters. The three-dimensional ﬁnite
diﬀerence or ﬁnite element calculations also have the capability to model these
low frequencies well, in those few regions where the velocity model is extremely
well known. Nonetheless, this eﬀect seems to indicate that special attention isneeded at large distances.
As a summary, we do see applications for synthetic seismograms generated
in ﬂat-layered velocity structures. This case study, based on a subset of sites in
Tohoku with Vs30>500 m/s shows both successes and opportunities for improvement. In this region with very complicated geology, the simple approximationof using a separate velocity structure appropriate for each site results in improvements in the standard deviations of the predictions compared to using asingle regional velocity model. Detailed site-speciﬁc models of the shallow geological structure brought mixed success and a modest reduction of the standard
6

−134−

deviation.

5

References

Anderson, J. G. (2015). The composite source model for broadband simulationsof strong ground motions. Seism. Res. Lett. 86, 68–74.
Anderson, J. G., R. J. Brune, J. N. Brune and G. P. Biasi (2017). Wave propagation and source models for synthetic seismograms compatible with strong
motion applications, Proceedings, 16th World Conference on Earthquake Engineering, Santiago, Chile, January 9-13, 2017, Paper No. 4023.
Brune, J. N. (1970). Tectonic Stress and the Spectra of Seismic Shear Waves
from Earthquakes, J. Geopphys. Res. 75, 4997-5009.
Ducellier, A., H. Kawase and S. Matsushima (2013). Validation of a newvelocity structure inversion method based on horizontal-to-vertical (H/V) spectral ratios of earthquake motions in the Tohoku area, Japan, Bull. Seism. Soc.
Am. 103 (2A), 958-970.
Goulet, C. A., N. A. Abrahamson, P. G. Somerville, and K. E. Wooddell(2015). The SCEC Broadband Platform validation exercise: methodology forcode validation in the context of seismic-hazard analyses, Seism. Res. Lett. 86,
17-26.
Kanamori, H. (1977). The energy release in great earthquakes, J. Geophys.
Res. 82, 2981-2987.
Kanamori, H. and D. L. Anderson (1975). Theoretical basis of some empirical relations in seismology, Bull. Seism. Soc. Am. 65, 1073–1095.
Kawase, H., F. Nagashima, K. Nakano and Y. Mori (2018). Direct evaluationof S-wave ampliﬁcation factors from microtremor H/V ratios: double empiricalcorrections to “Nakamura” method, Soil Dynamics and Earthquake Engineering(in press).
Kiran, K., S. Thingbaijam and P. M. Mai (2016). Evidence for truncatedexponential probability distribution of earthquake slip, Bull. Seism. Soc. Am.
106, 1802-1816.
Luco, J. E., and R. J. Apsel (1983). On the Green’s function for a layered
half-space, part I, Bull. Seismol. Soc. Am. 73, 909–929.
Rodriguez-Marek, A., G. A. Montalva, F. Cotton and F. Bonilla (2011).
Analysis of single-station standard deviation using the KiK-net data, Bull.
Seism. Soc. Am. 101, 1242-1258.
Somerville, P., K. Irikura, R. Graves, S. Sawada, D. Wald, N. Abrahamson,
Y. Iwasaki, T. Kagawa, N. Smith, A. Kowada (1999). Characterizing crustal
earthquake slip models for the prediction of strong ground motion, Seism. Res.
Lett. 70, 59-80.

7

−135−

6

Figures

EQ #3: Takahagi, MW=5.8, March 19, 2011
FKSH01
FKSH02
NIGH07

FKSH17

NIGH10
FKS007
FKSH09FKS025
FKSH05

NIGH12

FKSH08FKSH06°

FKS027

NIGH15
FKSH0737 N

FKS015

TCGH17

FKS014

GNMH09TCGH13

GNMH07
TCG010

IBRH14
IBRH16

TCGH14 TCG007
GNMH14
TCG011

IBRH18
IBR008
IBRH19SITH05°36 N

SITH09SIT007
SITH08SITH10
SIT006
SITH07
SIT012
YMN001
TKYH13
KNG011°139 E

140 ° E

Figure 1: Location map: event and stations.

8

−136−

141 ° E

EQ #3: Takahagi, MW=5.8, March 19, 2011
FKS015°37.0 N°36.9 N
FKS014°36.8 NTCGH13°36.7 NIBRH14

IBRH16°36.6 N°

140.2 E°

140.3 E

140.4° E

140.5° E

140.6° E

140.7° E

Figure 2: Detailed map of the epicentral region, showing nearby stations, epicenter of the main shock, immediate aftershocks, and the surface projection of
the fault plane used to model the earthquake. The fault is square, 7.96 km oneach side. The depth of the top edge is 4.0 km, and the hypocentral depth is
5.37 km. The fault strike is to the southeast, and it dips 48o to the southwest.

9

−137−

Figure 3: Cross-section view of aftershocks and the fault model. The view is
towards the southeast.

10

−138−

Figure 4: Earthquake hypocenter and aftershocks, projected onto the model
fault plane. View is downwards, from the southwest.

11

−139−

VM 1

vp

10
5
0
0

5

10

15

20

25

30

35

vs

5

VM 1
VM 2

0
0

5

10

15

20

25

30

35

0

5

10

15

20

25

30

35

0

5

10

15

20

25

30

35

Qp

1000
500
0

Qs

1000
500
0

Depth, km

Figure 5: Velocity models compared to 35 km.

12

−140−

VM 1

vp

6
4
2
0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

vs

4
2

VM 1
VM 2

0
0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

Qp

1000
500
0

Qs

400
200
0

Depth, km

Figure 6: Velocity models compared to 0.5 km.

13

−141−

Figure 7: IBRH16 observed and synthetic accelerations, calculated with VM4.

14

−142−

Figure 8: IBRH16 observed and synthetic velocities, calculated with VM4.

15

−143−

Figure 9: IBRH16 observed and synthetic displacements, calculated with VM4.

16

−144−

IBRH16, r flt = 20.6

Fourier Amplitude, cm/s

10

2

10 1

10 0

10 -1

Mean Horizontal Data
VM 4
VM 3
VM 2
VM 1

10 0

10 1

10 2

f, Hz
Figure 10: Observed Fourier amplitude spectra at station IBRH16, and average
predictions from each velocity model. These spectra have been smoothed, andshow the amplitudes of the vector sum of the horizontal components.

17

−145−

IBRH16, r flt = 20.6

PSA (5%), cm/s2

10

3

10 2

10

Mean Horizontal Data
VM 4
VM 3
VM 2
VM 1

1

10

-2

10

-1

10

0

10

1

Tn, s
Figure 11: Equivalent of Figure 10 for the pseudo-acceleration response of data
and average predictions from each velocity model.

18

−146−

Data/Model: Fourier Spectrum

10

1

10 0

10 -1
10 -1

10 0

10 1

Data/Model: PSA

f, Hz

10 0

10

-2

10

-1

10

0

10

1

T, s

Figure 12: Average misﬁt, by station, relative to synthetic predictions from VM
4. The color of the line depends of the distance, rf lt from the station to thenearest point on the fault. Colors follow the colors of the spectrum, with red
for the nearest station and violet for the farthest. The heavy black line is the
42 station average.

19

−147−

0

6

10

-3

FS(1 Hz)
PSA(1 s)

4

Slope

Slope

-0.005-0.01
FS(0.3 Hz)
PSA(3 s)20-2-0.015

-4
1

2

3

4

1

2

Case
5

3

4

Case

10 -3

10

10 -3
FS(10 Hz)
PSA(0.1 s)

Slope

Slope

5
0

0
FS(3 Hz)
PSA(0.3 s)

-5

-5
1

2

3

4

1

Case

2

3

4

Case

Figure 13: Residuals, by velocity model, as a function of distance at four diﬀerent frequencies.

20

−148−

ln[FS/Model]

1
Case 4
Case 3
Case 2
Case 1

0.5
0
-0.5
-1
-1
10

10

0

10

1

10

2

10

0

10

1

f, Hz

ln[PSA/Model]

0.5
0
-0.5
-1
-2
10

10

-1

Tn, s
Figure 14: Eﬀect of the velocity models on the mean residuals over 42 stations
as a function of frequency for Fourier spectra (top) and by period for response
spectra (bottom), by velocity model. The residuals for Fourier spectra areshown only to 12 Hz because we are not conﬁdent that the signal/noise ratio is
adequate at all 42 stations.

21

−149−

FS

1
Case 4
Case 3
Case 2
Case 1

0.5

0
10 -1

10 0

10 1

10 2

0

10

f, Hz

PSA

1

0.5

0
-2
10

10

-1

10

1

Tn, s
Figure 15: Eﬀect of the velocity models on the standard deviation of residuals
over 42 stations as a function of frequency for Fourier spectra (top) and by
period for PSA (bottom), by velocity model. For PSA, the solid and dashedlines show values of φ and φSS , respectively the intraevent and single-station
standard deviations measured at KiK-net stations, from Rodriguez-Marek et al.
(2011).

22

−150−

0.78

Mean Standard Deviation

0.77
0.76
0.75
0.74
0.73
0.72
0.71

Fourier Spectra
Response Sprectra

0.7
1

2

3

4

Velocity Model
Figure 16: Eﬀect of the velocity models on the mean standard deviation, theaverage over the spectrum of the standard deviations shown in Figure 15.

23−151−