Chapter 7

COULOMB STRESS CALCULATIONS

Here we calculate static stress changes caused by the displacement of a fault or dike or point source (we refer to these as ‘source faults’), derived from the strain field shown in previous sections. We resolve the shear and normal components of the stress change on grid node points or on specified ‘receiver’ fault planes. Receiver faults are planes with specified strike, dip, and rake, on which the stress changes caused by the source faults are resolved. The shear stress change (an increase or decrease) is dependent on the position, geometry, and slip of the source fault, and on the position and geometry of the receiver fault geometry (including its rake). The normal stress change (clamping or unclamping) alone is independent of the receiver fault rake.

We use the Coulomb failure criterion, Ds_f = D t_s + m’ D s_n, in which failure is hypothesized to be promoted when the Coulomb stress change is positive. Here, Ds_f is the change in failure stress on the receiver fault caused by slip on the source fault(s), D t_s is the change in shear stress (reckoned positive when sheared in the direction of fault slip), D s_n is change in normal stress (positive in unclamping of the fault), and m’ is effective coefficient of friction on the fault.

We’ll explore three kinds of receiver faults: (1) receiver faults listed in the input file with no slip, (2) focal mechanism files, and (3) faults optimally-oriented for failure. Their optimal orientations are a function of the regional stress, the stress change associated with the source fault, and the assumed friction coefficient.

On the following page are Figures 2 and 3 from King et al. [1994], which graphically present the Coulomb stress change resolved on vertical strike-slip faults parallel to the source fault (Fig. 2a); and resolved on optimally oriented planes (Fig. 2b) for a given regional uniaxial compression and friction coefficient. In Figure 3, the influence of the regional stress magnitude is seen on the orientation of the optimal planes, and on the stress change resolved on these planes. These figures were made in a primordial version of Coulomb.

Illustration of the Coulomb stress change (Fig. 2 from King et al [1994]). The panels show a map view of a vertical strike-slip fault embedded in an elastic halfspace, with imposed slip that tapers toward the fault ends. Stress changes are depicted by graded colors; green represents no change in stress. (A) Graphical presentation of equation 8 of King et al (1994), a “specified fault” calculation. (B) Graphical presentation of equation 13 of King et al (1994), for optimally-oriented strike-slip (“opt strike-slip”) faults.

Dependence of the Coulomb stress change on the regional stress magnitude, s^r, for a given earthquake stress drop, Dt (Fig. 3 from King et al. [1994]). If the earthquake relieves all of the regional stress (left panel), resulting optimum slip planes rotate near the fault. If the regional deviatoric stress is much larger than the earthquake stress drop (right panel), the orientations of the optimum slip planes are more limited, and regions of increased Coulomb stress diminish in size and become more isolated from the source fault. In this and subsequent plots, the maximum and minimum stress changes exceed the plotted color bar range (in other words, the scale is saturated).

7.1 Coulomb stress change on specified receiver faults (‘Specified faults’)

This is the simplest calculation, and is widely used by researchers. To resolve the stress, you need to specify the fault strike, dip, and rake of the receiver planes following the Aki & Richards convention, shown graphically in the next page of this manual.

1. Launch Matlab/Coulomb 3.1.

2. Choose Input > Open existing input file. Then in the Open input file window, choose the “input_files” sub-folder within the “Coulomb 3.1” folder, and select “Example-2(LL)-lonlat.inp”

3. Choose Functions > Stress > Coulomb stress change.

4. You will see the pop-up Stress control panel.

5. To calculate stress changes on specified receiver fault planes, click “Specified faults”. Coulomb averages the information on all input fault patches listed in the input file and puts the values in the boxes, but you can change them. Strike, dip, and rake are defined following the conventions of Aki & Richards (1980, 2002):

For practice, choose a strike/dip/rake of 360°/90°/0°, change the friction coefficient to 0.0, set the stress-change color saturation to ±5 bars, and hit ‘Calc. & View’.

The Coulomb stress changes default using our “Anatolia” color scheme. Other schemes (“Rainbow” or “Black & White”) can be changed by pulling down the Input > Preferences menu bar and clicking on “Color Map”. Then re-run the input file.

7.2 Using the strike/dip/rake/friction slider (Specified slip control panel)

To vary these parameters on the fly, hit the > button in the Stress control panel, and the Specified slip control panel pops up (see it above). Slide the balls to explore how the stress pattern changes.

7.3 Saving the graphic and numerical output of stress calculations

To save this graphic, File > Save As… > choose a .pdf format and rename it; its only 40 kb but is a full vector image. Stress change is calculated in the lower left corner of each grid square at the target depth specified in the input file. Calculated values of stress change may saturate (exceed the plotted range), so experiment with the slider to see more subtle features. The stress-scale legend is to the right of the plot.

Every time you click “Cal. & View” in the Stress control panel, a numerical output file called “dcff.cou” will be created or updated and saved in a sub-folder you designate in the Preferences menu. Remember to rename “dcff.cou” if you want to save the numerical file; otherwise they will be overwritten.

7.4 The importance of the receiver fault geometry in Coulomb modeling

In the figure below, from Lin & Stein (JGR, 2004), the source is an idealized M~7.9 right-lateral southern San Andreas rupture. The most familiar stress pattern is the case of source and receiver faults with the same strike, dip and rake (a). But look how strongly the pattern changes when the receiver faults are nearby thrusts (b-c) or left-lateral faults (d). Thus, a stress shadow for one receiver may be a stress trigger zone for another receiver fault. Note also that we tend to assume a high coefficient of friction (~0.8) for continental thrust faults, moderate friction (0.4) for strike-slip or unknown faults, and very low friction (>0.2) for major transforms, such as the San Andreas.

7.5 Adding Coulomb stress to an overlay plot and viewing it in 3D

1. Launch ‘Example-SFBayArea.inp’. Overlay > Coastlines. Open the ‘coastline data’ folder in the Coulomb30 folder, and select ‘california_coastline_di_neg.dat’ Answer the pop-up question, ‘no’ (some datasets treat western longitudes as negative; others do not). Now, Overlay > Active faults. Choose ‘California faults_longlat_dat’.

2. You can now add any kind of stress changes to this plot without losing any of the overlays. Just hit Functions > Stress > Coulomb stress change, and select a stress component and option, and you will see a map view image such as that below left.

3. After completing any Coulomb stress calculation made over the grid, such as Coulomb stress change, Shear stress change, or Normal stress change, just return to the MATLAB window and type in "coulomb_3d_view". This calls a plugin (which can be found in the plugin folder in Coulomb30 folder). It executes a 3D view that you can manipulate and save as a pdf file, such as that above right.

7.6 Coulomb stress cross-sections

1. Choose Input > Open existing input file > Example-2(TH).inp. Next, choose Functions > Stress > Coulomb stress change. Change the color saturation to 3.0 bars, and hit ‘Calc. & View’ in the Stress control panel. After the mapview plot appears, hit the Cross section button on the Stress control panel, which displays the cross section parameters from your input file (left plot below). Click “Calc. & View” and you will see the image below left. Now change the saturation to ±10 bars and chose interpolated shading, and hit “Calc. & View” again (right panel below).

Mosaic, ±3 bars Interpolated, ±10 bars

2. When building a cross-section, the “Start (A)” point should always be on the left side of the “Finish (B)” point (see below right). You can overlay seismicity in cross-sections if you first plot the mapview with earthquakes and then make the section. The default setting plots seismicity within ±20 km of the section line. To change this, print EQPICK_WIDTH=5 for ±5 km, etc, in the command window.

2. You can also change the dip angle of the cross section by entering a new value in the “Dip” box in the above panel. The blue dashed line is the depth at which stress was sampled in the map view. The red line is the intersection of the fault plane with the cross-section line. These can be seen in map view, which has a blue cross-section line, a red fault perimeter, a green line where the fault projects to the ground surface, and a black line where the fault intersects the depth at which stress is being sampled. You can change the depth, and many other parameters, in the Stress control panel.

4. You can calculate the maximum or mean values of Coulomb stress changes between various depth ranges. Click “Depth range” and then “Cal. & View” in the Stress control panel, you will see Depth control panel. In the Depth control panel, enter depths of the top and bottom surfaces, as well as the depth increment, for which you want to perform the comparisons. You can calculate either “Maximum values” of Coulomb stress changes over the given depth range (below left), or the “Mean values” over the range (below, right). The numerical output file shows the maximum value at each grid point. The maps below were made using the left-lateral strike-slip fault input file, Example-2(LL)-lonlat.inp.

Max stress change over a depth range Mean stress change over a depth range

5. You can plot the orientation of the ‘specified fault’ strike in map view. Choose Input > Open existing input file > Example-2(LL).inp. Next, choose Functions > Stress > Coulomb stress change. Check the “Strike line” box in the Stress control panel, and for practice, change friction to 0.0, shading to interpolatedand then “Calc. & View”, resulting in the following image on the left. Notice that the lines strike 41° as specified. Stresses and slip lines are plotted on the lower left corner of each grid box.

7.7 To display the principal stress axes rotated by the earthquake stress change

Coulomb makes the tensor addition of the earthquake stress change, which diminishes with distance from the source fault, and the regional stress, which is assumed uniform, to calculate the total principal stress axes. Once the total principal stress axes are determined, Coulomb uses this and the assumed friction coefficient to determine the optimum planes.

Choose Input > Open existing input file > Example-2(LL).inp. Then, Functions > Change parameters > Grid, and change the grid to 5.0 x 5.0 km so the axes are not too crowded (this is not a requirement; you can choose any grid spacing of interest). Now, Functions > Stress > Coulomb stress change. In the pop-up Stress control panel, choose “Opt. Strike S.”, “Principal Stress”, and “Interpolating”, and click “Calc. & View”. You will see the image on left below. Using Tools > Zoom In, you will see the right image below.

Normal view with principal axes Zoomed view

The axes are not scaled by their stress magnitudes, since too many would be invisible. Notice that the axes rotate in 3D, not just in the horizontal plane, and so s₃ deviates from being vertical close to the source, as a result of the stress changes imparted by the earthquake. This means that s₃ can be seen in the lower left corner of the zoomed image.

7.8 Using the regional or ‘tectonic’ stress in optimally-oriented stress calculations

The regional stress is used only when you choose Functions > Stress > Coulomb stress changes AND choose one of the optimally-oriented (“Opt.”) stress changes in the Stress control panel. Otherwise the regional stresses are ignored (see below). See King et al (BSSA, 1994) for more on this topic.

1. All of the example input files contain simple regional stress fields. The sign convention for S1, S2, and S3 is positive in compression [bars]. Although you can assign axes that are neither horizontal nor vertical (in other words, they plunge at angles other than 0° and 90°), and you can assign a vertical gradient in the stress axes, such that the stress magnitude increases (or decreases) with depth, we advise you to keep it simple. We almost always make the greatest principal stress, S1, horizontal and about 100 bars, and the intermediate principal stress about 0 bars.

2. A reasonable regional stress must promote slip on all the faults in your area. You can use the P and T axes from focal mechanisms of earthquakes off the main faults, hydrofrac or borehole breakout information, and fault striae inversions to assess the strike of the greatest and least principal axes. The World Stress Map (http://www-wsm.physik.uni-karlsruhe.de/) may also be useful for this purpose. But most important, ask yourself if the stress tensor you are proposing is consistent with the fault types and orientations in the area of study.

3. You can build principal stress axes (which must be mutually perpendicular) by selecting Functions > Tools > Principal axis calculator. You can also visualize the principal axes by Functions > Stress > Coulomb stress change, and then by checking ‘Principal axes’ in the pop-up Stress control panel.

7.9 To display Coulomb stress changes on the optimally oriented strike-slip, thrust, and normal faults (with the fault dip fixed) in map view

“Optimally-oriented” planes are those oriented such that the Coulomb stress change resolved on them is most positive (in our plots, most red). Here one assumes that aftershocks have a range of orientations, and will be most abundant at the depth, location and with the geometry that most promotes failure. So at every point on our grid, we find the receiver plane along which the resolved Coulomb stress change is most positive. If the regional stress were zero, there would be a positive Coulomb stress change on some plane everywhere. But at large distances from a source fault, the earthquake stress change will be much smaller than the regional stress, so the planes rotate very little.

1. Choose Input > Open existing input file> “Example-2(LL).inp”. To calculate optimally oriented faults, one must specify the regional (sometimes called, tectonic) stress field in the input file. “Example-2(LL).inp” contains a sample 3D regional stress field. If you use your own input file, make sure the regional stress tensor is properly defined. To make sure that the directions of your regional stress axes are geometrically compatible (i.e., orthogonal to each other), choose Functions > Tools > Principal axes calculator. You will see a pop-up Coordinate axes for principal stress window (below).

To check the compatibility of regional stress axes, first enter the azimuth (strike) and plunge (dip) of the greatest principal axis, S1 azimuth and S1 plunge. Click “Calc”. Then enter the plunge of the intermediate principal stress axis (“S2 plunge”), which should stay within the shown permissible range of values (“Plunge range”) to ensure that three principal axes are mutually perpendicular, which is obligatory. Click “Calc” again, and you will see four sets of combinations for the three principal axes, all of which are consistent with the entered values.

2. Choose Functions > Stress > Coulomb stress change. You will see the following in the pop-up Stress control panel:

In Stress control panel, the alternatives are “Opt Thrusts”, “Opt Normal”, and “Opt Faults”. ‘Opt fault’ is an exceedingly slow grid search over the entire focal sphere for the receiver planes at every grid point, and so you need patience in waiting for the completing of the calculation. Here, choose “Opt Strike S.” and click on Strike lines; then click “Calc. & View” and you will see this in the interior of the plot:

1. Explanation of the Slip lines. For the “Opt Strike S.” option, the receiver faults dip vertically and their strikes are a function of S1 (a₁) and the friction coefficient m (FRIC in Coulomb). The right-lateral and left-lateral optimum planes are plotted at every grid point. Their internal angle is a function only of the friction coefficient; they are orthogonal for zero friction and acute for high friction. You can see large rotations of these planes very close to the source fault, because the earthquake stress change is about as large as a the assumed regional stress, ~100 bars. Far from the source, the slip lines are parallel and their orientation is a function of the regional stress in the input file, because at these distances from the source, the earthquake stress changes are negligible. For the thrust and normal faults, their dip is a function of the friction coefficient, and their strike is perpendicular to S1 (a₁) and S3 (a₃), respectively. The relationship of the angle, b between S1 and the dip (for strike-slip faults) or strike (for dip-slip faults) is tan 2b = 1/FRIC. See King et al. [1994], which can be accessed online from the Coulomb web site at

http://quake.wr.usgs.gov/research/deformation/modeling/papers/landers.html

4. For a friction coefficient of zero, the optimum planes would be orthogonal and we could compare them to aftershock focal mechanisms unambiguously. Try checking this by selecting Functions > Change parameters > Coeff. Friction. In the pop-up Coeff. Friction window, type 0.0 and click “OK”. Then, Functions > Change parameters > Grid, and change the grid to 5.0 x 5.0. (For normal and thrust planes, the strike for the two planes is the same and their dip is a function of the friction coefficient). Then hit Functions > Stress > Coulomb stress change. In the Stress control panel, check ‘Opt. Strike S.’, and ‘Strike line’. An excerpt is shown below left. Repeat the process and in the Stress control panel, change the friction to 0.8 and hit the tab. You’ll see the image on the right.

Friction coefficient = 0.0 Friction coefficient = 0.8

The stress changes are larger and more asymmetrical about the fault plane for high friction, and the rotation of the optimal planes near the source fault are a little greater. These aspects are discussed in King et al [BSSA, 1994].

5. Strike lines for optimally-oriented thrust faults. Choose Input > Open an existing file and choose ‘Example-2(TH).inp’. Then choose Functions > Stress > Coulomb stress change. In the Stress control panel, check ‘Opt. thrusts’, ‘Strike line’, and ‘Tiled’ image shading, and hit ‘Calc. & View’. You will see an image like this:

Here the two strike lines are parallel to the mapview (they would not be so in cross-section), so only one is visible. The rotations are greatest in the deep red and blue areas, where the stress change is large relative to the assumed regional stress. In this example, the fault is not quite orthogonal to the assumed regional stress.

6. In addition to “Opt Strike S.”, Opt Thrust”, and “Opt Normals”, you can also calculate optimally-oriented receiver faults on which Coulomb stress is maximum while not restricting fault types. Click “Opt Faults” and “Calc. & View” in the Stress control panel. You will see the following pop-up window:

You are requested to enter interval of grid research (°) in the Grid interval window. Small intervals take time to calculate. We recommend an interval of 10° or greater, at least until you’ve settled on all the parameters you want to use.

7.10 To resolve stress changes on receiver faults (‘Calc. stress on faults’)

This is perhaps the simplest Coulomb calculation, and is widely used in research. The regional stress information is not used, regardless of whether it is included in the input file. Coulomb permits one to show the resolved stresses graphically and output the numerical results.

Open ‘Example-1.inp’ in an editor; notice that there is a source fault with slip, and two receiver faults (called receivers because they lack slip). These two faults are automatically broken up in to 5 patches by placing the number ‘5’ in the ‘#’ column. This permits a stress calculation at the center of each patch on a plane. This is much easier and faster than inputting 10 receiver faults.

Now, Input > SFBay_rake.inp; Functions > Stress > Calc. stress on faults. A Stress on faults control panel will open to let you choose the stress component that you would like to resolve on the faults. A specified uniform rake of 135° was used in the plot below, with the (blue) stress drop on the source fault and one of the receiver faults, and the (red) stress increase on the other receiver patches, which varies along strike.

If your input file used the (rake, net slip) format rather than the (rt.lat. slip, reverse slip) format, then you could use the ‘Coulomb for individual rake’ option. To switch how the input file is interpreted, it searches for the word, “rake”, so don’t misspell it. In the input folder, the full 1400-patch California Reference Database is included in various formats. ‘CA-fault-database-rake.inp’ (and ‘CA-fault-database-rake.mat’) are formatted by rake for this purpose. Because it is a large file, it takes time to calculate. Here are the first two sources in the database input file, with the right-lateral (rake rake=180° ) Green Valley reverse (rake=90°) Mount Diablo thrust faults:

You can also use a subset and plot it, since all the faults have their associated names. Remember to make sure the ‘#fixed=’ in the third line of the input file is changed to reflect the number of faults in the input file (‘#fixed=faults+receivers).

A more complex case is represented for the 1995 Kobe earthquake. Both a uniform-slip and variable slip source model input files are included in the input file folder, in addition to about 50 surrounding faults with various orientations and dips. Input > Preferences > Cartesian; then Input > ‘Kobe-uniform-slip.inp’; then Functions > Stress > Calc. stress on faults. Note that a numerical output file is created each time you make a Calc. stress on faults run, which is very useful for further analysis. Two stress components are shown below. Experiment with different viewing angles and azimuths using the numerical input window; using the rotation tool with files this large is difficult due to the time lag.

Stress from the 1995 Kobe variable-slip model resolved on surrounding faults

Notice that on the Kobe rupture itself, there are both red and blue patches. The Coulomb stress change is not uniformly negative (blue) because of the irregular coseismic slip which leaves spikes of stress increase, the presumed sites of along-fault aftershocks.

7.11 Calculating stress change on a subdivided (tiled) receiver fault

Often one would like to determine the distribution of Coulomb stress change on a large receiver fault. This can be easily done and visualized in Coulomb 3.1.

1. Open ‘Example-1.inp’ in a text editor and change the ‘5’ in the # column for ‘receiver-2’ to 1, so it's now a single patch. Then launch ‘Example-1.inp’ in Coulomb. Functions > Tools > Taper or subdivide fault slip, type fault 3 in the fault # box, and split the source into 8 x 4 patches, and hit OK (below left). So now you have subdivided it into 32 patches.

2. Functions > Stress > Calculate stress on faults. You can zoom in and spin the view, as shown above right. The Coulomb stress change is colored and the rake vectors are plotted in blue. You can change the rake in the Stress on faults control panel.