Chapter 8

ADVANCED STRESS FEATURES AND MAGMATIC SOURCES

At this point we assume that you are now familiar with the typical steps and menu options in Coulomb, and our instructions will be more abbreviated, but still accurate.

8.1 Coulomb stress changes resolved on earthquake nodal planes

1. If there is focal mechanism information available, you may want to test whether the fault plane associated with each earthquake was brought closer to Coulomb failure. This is presumably a more stringent test than determining if earthquakes occurred in the red regions, although typically only the fraction of the catalog contains reliable focal mechanisms.

If the friction were 0, this would be a straight-forward test because the shear stress change on the orthogonal nodal planes is identical (eg., Parsons, JGR, 2002). But for non-zero friction, the optimum Coulomb planes are not orthogonal, and so they will be different. In fact, one plane may have a stress increase and the other a decrease. You can sometimes identify the fault plane based on independent information, and there are published approaches others have taken, such as comparing the percentage of stress increases on seismicity in the region before and after a mainshock (Hardebeck et al, JGR, 1998; Ma et al, JGR, 2005), or by using seismicity alignments (Seeber & Armbruster, Nature, 2000).

2. Using the same ‘Example-SFBayArea.inp’ you used in the last section, Overlay > Clear overlay data from memory > Clear earthquake data (because the file used did not contain focal mechanisms). Now, Overlay > Earthquakes > In the Earthquake catalog formats window, select ‘NCSN fpfit (focal mech)’, and then choose the dataset with this name from the folder. Earthquake circles plot on the map.

3. Functions > Stress > Calc. stress on nodal planes. If the focal mechanism file only includes one set of nodal planes, Coulomb asks you if you want it to generate the other set, so answer ‘yes’. A window appears to ask you which set (Nodal set 1 or set 2, which are usually arbitrarily ordered in a catalog, or Shuffle—random selection of set 1 or 2). It will assign colors to each earthquake according to how much the nodal planes were brought closer to or farther from failure. You should see this plot:

You can plot each set (Nodal plane 1 or Nodal plane 2) separately, shuffle them randomly, or you can edit the catalog in an editor to make set 1 the likely fault planes. The numerical output file also calculates the maximum (most positive) Coulomb stress on each plane, and gives the corresponding rake for the maximum.

8.2 Coulomb stress at a point, from lat/lon, x/y, or a mouse click

Sometimes you want to probe the calculated stress changes at a given point, rather than over a regular grid or on receiver faults. This is easy to do in Coulomb by using the mouse or by entering a lat/lon or UTM coordinates.

1. Functions > Stress > Calc. stress at a point. A Stress calculation on a point window pops up. You can enter the coordinates and depth of interest, and you can also enter the receiver fault geometry. The output is numerical, and appears in the Output window, with blue numbers when negative and red when positive (see below).

2. When using the Calc. stress at a point option, you can change the assumed friction coeffient in the same manner that you can for other stress calculations, by choosing Functions > Change parameters > Change Coefficient fric., and entering the new number. If you want to use x/y values, change to Cartesian coordinates in the Preferences menu. Then return to the Calc. stress at a point selection and use the Stress calculation on a point window. The lat/lon vs. Cartesian options are highlighted above in the red box.

8.3 Point shear sources for small quakes or curved faults (Kode 400)

1. If you specify Kode 400 rather than 100 in the input file, Coulomb will interpret it as a point source with rt.-lateral and reverse slip components. The location of the point source is defined by the rectangle in the input file. The point is located at the center of the rectangle: the midpoint of the start and end positions, and at the mid-depth of the top and bottom positions. The plane defined by this geometry is used for strike and dip of the point source. This is a simple way to represent small shocks.

2. Another use of Kode 400 sources is if you have a very crenulated or curved fault, the overlap and underlap of the rectangles can be avoided by making a mesh of point sources (shown schematically below left). But this is demanding, and don’t sample the stress closer to the fault than the point spacing, or you will see numerical noise.

3. We provide an input file with a single point source, ‘Example-Kode400.inp’. Note that the rt.-lateral and reverse units are in potency (m³), not in slip (m), since a point has no area (see above right, from open & edit input file). Potency, also called the ‘geometric moment,’ is the seismic moment/G, where G is the shear modulus (see boxed equation in Chapter 3 to calculate G from E and PR).

8.4 Crack opening or closing (Tensile Kode 200 and 300 sources)

Many magmatic intrusion problems involve a tabular body, such as dike or sill, undergoing expansion or contraction. Here we introduce the deformation and stress change caused by a tensile component of displacement, such as the opening or closing of a tabular or crack-like body, using the equations of Okada [1992]. For simplicity, the input file format for the Coulomb has only two slip components, labelled ‘rt.lat’ and ‘reverse’. So, when you change the Kode of the source to 200 or 300, the first of these slip components takes on a tensile-displacement interpretation by the program, as described below.

1. Kode 200 permits tensile (opening) slip and strike-slip to be input. Put lateral slip in the ‘rt.lat.’ column and tensile slip in the ‘reverse’ column (dike opening is positive; unit = meters), ‘Example-6(200Kode).inp’. Here’s how a Kode 200 source appears in the Functions > Open & edit input file window,

2. Kode 300 permits tensile (opening) slip and dip-slip to be input. For Kode 300, put tensile slip in the ‘rt.lat.’ column (again, dike opening is positive), and put reverse slip in the ‘reverse’ column.

3. If the element has three components of slip (strike-slip, dip-slip, and tensile-slip), you can accomplish this by including both Kode 200 and a Kode 300 sources in the input file as shown in the file below, ‘Example-7(200-300kode).inp’. To the right of the ‘bot’ values you can add comments that help you keep track of the source or receiver name; this is true for all .inp files.

1. Define the dike geometry (endpoints, dip, and top/bottom depth) as you do for faults. Unfortunately, the way we define sources and receivers does not permit one to put in a strictly horizontal plane, since the top and bottom depths would be identical, but you can put a tiny dip (~1°), and it will work fine.

5. Uniform opening or closing of a dike is likely to be unrealistic. Just as you can for faults, you can taper the displacement toward the edges by Functions > Tools > Taper.

6. Select the file, Example-8(dike).inp, and then choose Functions > Displacements > 3D Vectors. You will get the image shown on the left. Now, Functions > Change Parameters > Vertical Exaggeration > 25000, and change Calculation depth to the surface (0.0 km). You will get the image on the right. Experiment with view angles.

3D Vectors at 10 km depth for a dike 3D Vectors at the free surface

8.5 Point source of deflation or inflation (Kode 500)

We have implemented the formula of Okada (BSSA, 1992) for a buried point source of expansion or contraction. Of course, the calculated stress changes very close to such a point source are infinite and are thus unreliable, but at larger distances, this source is useful for magmatic intrusion problems, such as studying how an intrusion event may bring other faults or dikes closer to failure. An example of this is Nostro et al., Two-way coupling between Vesuvius eruptions and southern Apennine earthquakes by elastic stress transfer, JGR, 103, 24,487-24,504, 1998, which is on our web site at http://quake.usgs.gov/research/deformation/modeling/papers/nostro.html. As with the dike or sill source, we have created two kodes for mixed volumetric and shear displacements.

1. The unit for deflation or inflation is potency (m³), the volume change associated with the evacuation or expansion of a spherical magma chamber. Expansion is positive. Potency, also called the ‘geometric moment,’ is equivalent to the seismic moment/G, where G is the shear modulus. As discussed previously (see boxed equation in Chapter 3), in Coulomb one inputs Young’s modulus (E) and Poisson’s ratio (PR), from which G can be computed. The units for the strike or dip slip is meters, as usual. One designates the source as Kode 500.

2. Kode 500 permits both a dike and point source of expansion to be input. Put the tensile opening in the first (‘rt-lat.’) column (units are m³; right-lateral is positive), and put the point source in the ‘reverse’ column (units are m³; expansion is positive).

3. The location of the point source is defined by a plane in the input file. The point is located at the midpoint of the start and end positions, and at the midpoint of the top and bottom positions. The plane defined by this geometry is used for the dike opening. So if you only calculate the point-source expansion, the length and width of the plane are unimportant; only if you also want to consider opening or closing of the planar dike do these dimensions get used in the calculation. The point source is indicated graphic output files as a small black circle.

The following figure was made using ‘Example-9(Kode500).inp’ with two different Displacement options selected. The surface projection of the dike is the short green line.

Vertical displ. (color & contours) 3D image drape

In the example below, we have chosen Stress > Coulomb stress change > Opt Normal with a 1 x 1 km grid, friction=0.8, depth of 8 km, and color saturation of ±0.25 bars:

Coulomb stress change (bars) on optimally oriented normal faults

8.6 Displaying Coulomb Displacement, Stress and Overlay graphics in Google Earth

Google Earth offers wonderful opportunities to drape Coulomb files on satellite images and DEMs. You need to download and install Google Earth on your computer before following these steps. In addition to this, you will need Adobe Illustrator or any reasonable graphic application.

1. Make sure that you have selected the lon. & lat. Coordinate option in the Preferences menu.

2. Run any horizontal displacement function, such as Horizontal displ. (vectors) or Horizontal displ. (wireframe). Add the Overlay options that you would like in the plot. When the final plot appears in the Coulomb 3.1 window, type “coulomb2googleearth” in the MATLAB command window, and hit return. You will get messages in the command window once files are saved. Then, you can drag the kml files onto the Google Earth icon on your computer dock or desktop, or you can open the files from the Google Earth menu.

This file was made using the input file, ‘Landers-uniform-slip.inp’ and Displacement > Horizontal Displ. (Vectors).

3. For color gradient images such as are created in the Coulomb stress change, Shear stress change, Normal stress change options, or the Vertical displacement (color) options File > Save As... pdf using the commands, renaming them.

4. Open the .pdf files in Illustrator or the graphics software you are using, and remove all the text and lines, and keep only the color (raster) image. Save the color image as .jpeg file, and name it ''coulombmap.jpeg'' in the Coulomb directory.

5. Type ''coulomb2googleearth'' in the command window. This calls a script that is run by Coulomb; this is not a feature of the normal pull-down menus. You will get the overlay kml file.

In the plot above, a Coulomb stress change image was added. The DEM exaggeration (3:1) was maximized in Google Earth with an oblique viewing angle from the south. Google Earth also lets you control the translucency of the jpeg.

6. The full plugin script can be found in the plugins folder in the Coulomb30 folder:

8.7 Digitize Polygon Function (to select and create a lat/lon file of any area, such as a Coulomb stress-change or deformation contour)

The Digitize Polygon function is useful for comparing a zone of Coulomb stress change or deformation with aftershocks or the seismicity rate change as calculated from the independent Matlab program, ZMAP.

1. Make sure that you have selected the lon. & lat. Coordinate option in the Preferences menu.

2. Run any map view function (not 3D) in Coulomb, and then type “digitize_polygon” in the MATLAB command window, and hit return. In the example below, we have used the input file, ‘Example-SFBayArea.inp’ and selected a Coulomb mapview. Below we have digitized the -3 bar contour.

3. The mouse becomes a digitizing cursor in the Coulomb main window. You will see a large “+” over the map (above right). To digitize, make multiple left-button clicks to circumscribe an area. Finalize and close the area by using a right-button click; the graphic will then show the completed area (above right). A window (left below) asks if you intend to digitize multiple areas; hit ‘no’ if this is it. The digitized file is saved in the "coulomb" folder. The first 8 lines of the file, ‘polygons1.dat’ is shown below at right; notice it is in lat/lon coord’s with the convention that western longitudes are negative.

4. In ZMAP, you can move the file to the ZMAP folder and then read it from the area selection function.

5. The ZMAP software, by Prof. Stefan Wiemer at ETH Zurich, is available on http://www.earthquake.ethz.ch/software/zmap