It is extremely useful to gain an even rudimentary idea of the attenuation in earthquake intensity with distance from the focus. Although it is not possible to get a complete picture of MSK intensity using simple computational methods, such methods can provide a heuristic idea of attenuation. The main goal of this lesson is to learn about seismic intensity attenuation models. You will use Perl scripts for calculating seismic attenuation and GMT to plot estimates of maximum earthquake intensity for an area and period of time of your choice. A secondary goal of this lesson is to become familiar with more complex Perl scripts for data processing and GMT scripts that contour data on a geographic map.
Numerous methods and formulae have been proposed for estimating earthquake intensity and related measured of motion (peak ground acceleration, peak velocity, root-mean-square ground acceleration). A brief summary of sum of these methods can be found in the new Manual of Seismological Observatory Practice:
Estimating seismic intensities is a complex undertaking. Here we will use one of the simplest, and as a result among the least accurate(!), approaches to estimating seismic intensity.
Methods for estimating the attenuation of MSK intensity with distance from the focus were first proposed by Kövesligethy (1906), and often revised (e.g., Blake). A typical model of attenuation of intensity looks like:
where magnitude is the surface wave magnitude, focal_distance is the great circle distance between the earthquake focus and the point of interest (taking into account the depth of the focus). The constants (1.5, 3.0) vary in many formulations and may be regionally specific to some extent. The equation is grossly well behaved. Intensity will increase with earthquake magnitude and decrease exponentially with distance from the focus. Note that this formulation is quite similar in form to some descriptions of peak ground acceleration (e.g., Joyner and Boore):
Another thing to note about the formula used for attenuation of intensity is that it is radially symmetric. Of course, this is a serious assumption: most actual intensity maps are not radially symmetric.
In application, one is immediately faced with the problem of calculating the great circle distance. This is the angular distance between the earthquake epicenter and the point where intensity is estimated. Why not just use the distance formula? Because we live on a roughly spherical earth. The straight line distance between points on a sphere is through the sphere itself. We are interested in the distance between two points on the surface of the earth. Here is a perl subroutine for calculating the distance between two points on a sphere.
Of course this leads to some errors since the earth is not a perfect sphere. The equation for the change in earth radius as a function of latitude is:
where the radius of the earth at the equator is 6378.139 km, f=1/298.247, and phi is latitude. It is possible to reduce the error in the the perl subroutine (about 2%) by using a latitude between the two points of interest, or going with more sophisticated approximations of the figure of the earth.
With a little work, and an earthquake catalog, the maximum intensities for a site can be estimated for some set of earthquakes.
The steps involved in producing the above map were:
1) Select an area of interest
2) Find large magnitude earthquakes in the catalog for this region of interest:
CNSS worldwide earthquake catalog
Southern California Earthquake Data Center
3) modify the perl script: msk_intensity.pl to calculate intensities for your region
4) modify the gmt script script8.gmt to plot the MSK intensities for your region.