I Created an Earthquake

SPECFEM3D Globe: SimCity for Adults

Benjamin Saadia
7 min readJan 4, 2022

It’s 9:10 PM on July 22, 2020. Eighty-eight Alaskans in the small village of Perryville are preparing for bed. Just two minutes later, a powerful 7.8 magnitude (M7.8) earthquake strikes 100 kilometres south of the Alaskan peninsula (see Figure 1 below). Strong shaking with a peak ground velocity of 16.83 centimetres per second is observed by a nearby seismic station. The next morning, locals report minor damages to the docks and the road leading to the harbour.

Figure 1: Location of the M7.8 July 22nd, 2020 earthquake epicentre (yellow star) near Perryville (red diamond). It originated due to thrust faulting near the Alaska-Aleutian subduction zone (solid white line), where the Pacific and North American Plates converge at a rate of 64 millimetres per year. Map modified from the USGS.

Preparing for Earthquakes

In this case, the proximity of the earthquake to human settlements meant that the damage was minimal. However, the impacts of earthquakes for those living in earthquake-prone regions can be enormous. Pacific-North American plate convergence has produced several historic high-magnitude seismic events, including the M7.9 Denali earthquake in Central Alaska, which resulted in ground offsets as high as 9 metres and an estimated $20 million U.S. dollars damage to infrastructure (United States Geological Survey, 2020).

While we cannot predict exactly when and where earthquakes will occur, we can use several streams of evidence to prepare for the impact these will have on our communities. I am a geoscientist with an engineering background, so this article presents a study I conducted that contributes to our understanding of seismic hazard and the methods needed to evaluate it. Seismic hazard is a complex topic that lies at the intersection of geoscience, engineering, urban planning, politics, and sociology. This is to say that the methods presented in this article are but one line of evidence that should be considered in our quest to understand how to prepare our communities for earthquakes.

Creating Earthquakes

I will admit that the title of this article is provocative if not misleading. Scientists are not supervillains, neither in desire nor capability. But in the age of abundant computing power, you are able to feed your God-complex in a far more accurate, and useful way than SimCity’s disasters allowed.

SPECFEM3D Globe (SF3D) is an open source software that can simulate seismic wave propagation on a global and regional scale (Komatitsch et al., 2015). While it’s not as fun as unleashing an earthquake on your SimCity creation, it is completely free and open-source. Given a set of inputs, it can simulate the seismic response at any point on Earth including the effects of topography and bathymetry, composition and density of the Earth’s crust, and ellipticity, rotation, and gravity.

This is powerful. If we can accurately simulate the seismic response at any point on the planet, we can predict which communities might be at risk if an earthquake occurs at a given location with a given mechanism. Note that this is not the same as predicting when and where an earthquake will occur. We are simply predicting what a seismometer would record at a location, given a particular earthquake.

Simulating the July 22nd, 2020 Perryville Earthquake

In this study, I simulated the seismic response to the July 22nd, 2020 earthquake from all of the stations that are a part of the Canadian National Seismic Network (CNSN) (Geological Survey Of Canada, 1989). The purpose of the study was to develop a generalized methodology for simulating seismic wave propagation at a regional-scale in Canada. The M7.8 Alaska earthquake is an excellent test-case given its high magnitude and proximity to Canada.

The precise technical information is as follows. The magnitude 7.8 Mw earthquake occurred at 06:12:44 UTC with its epicentre located at 55.072°N, 158.596°W (±5.4 km) and at a depth of 28.0 km (±1.9 km). For more technical information, see this summary from the USGS (2020).

Methods

SF3D 7.0 was used for generating synthetic seismic signals which were later processed using Obspy, an open source seismic analysis Python package (Beyreuther et al., 2010). Several processing steps were required for generating a STATIONS file containing the CNSN station names, latitude, longitude, and elevation in the correct format for SF3D. This file tells SF3D where on Earth to calculate the results. Because I wanted to compare the real measurements at the time of the earthquake to the model, I chose to use a list of locations that represent real seismic stations. The moment tensor solution provided here by the USGS tells SF3D what the earthquake source is like, physically and mechanically.

The remaining processing steps were carried out on both the real and synthetic seismogram data for the July 22nd earthquake. The real seismic data was initially processed by removing the instrument response from each station. This is necessary because each real recording was gathered using an instrument with a known sensitivity to certain frequencies. SF3D assumes a perfect instrument, so I needed to remove these known effects from the real recordings to make a good comparison.

Both datasets were then filtered using a bandpass filter with a minimum and maximum frequency of 0.023 Hz and 0.025 Hz, respectively. These limits were chosen given that the earthquake source’s half-duration period is known to be 41.08 seconds, or 0.0243 Hz. This does two things. Firstly, it filters out any numerical artifacts in the simulated seismograms. Secondly, it removes the effects of both random instrument noise (whatever is left after removing the systematic noise of the instrument response) and other seismic activity. For example, a seismometer located in an urban area may pick up vibrations caused by nearby traffic in the 20 to 50 Hz range. By removing these higher frequency waves, only the waves caused by the earthquake (and other sources in the 0.023 to 0.025 Hz range) will be present in the data.

For a complete summary of the workflow, refer to Figure 2 below.

Figure 2: A flowchart showing the complete methodology used in the study. This is generalized so that it can be used for any study involving the CNSN database with minor modifications. Note that the CAC cluster refers to a supercomputer I had access to for the duration of the project.

Results and Discussion

In the study, I generated results for the entire CNSN database. However, given that I am located in Kingston, Ontario, I figured that I would present the results at KGNO, the seismic station located right near my office. See Figure 3 below for these results.

Figure 3: Comparison of seismograms at station KGNO in Kingston, Ontario. The top row represents the unfiltered data, while the bottom row is filtered between 0.023 and 0.025 Hz. Red traces are the real, measured seismograms and blue are the synthetic seismograms generated in SPECFEM 3D Globe. Each column is a different measurement direction, Z being up-down, E being East-West, and N being North-South.

From Figure 3, the importance of filtering the data around the source frequency is clear. Visually, the unfiltered results (top row) do not match up. My assumption is that this is due to the effects described earlier: numerical artifacts in the simulation and ambient noise in the real measurements. Once filtered, the power of SF3D becomes clear.

Once filtered, the signal polarity and phase are very close. The amplitudes are not identical however, which is probably a result of s20_rts, the global velocity model that was used in the simulation. The velocity model defines the makeup of the earth material through which the waves propagate. The better the velocity model, the better the results should be. This model is global, meaning it does not take into account local variations in rock properties. While the magnitude of displacements are not equal, they are synchronized. For example, in the filtered Z-component plot (bottom left) in Figure 3, the amplitudes of the synthetic (blue) and measured (red) signals increase to a maximum at about 1100 seconds. The same phenomenon can be observed in the E- and N-component filtered plots (bottom middle, and bottom right).

See Figure 4 below for an animation of the July 22nd earthquake, along with the simulated signal at KGNO. Turn on your volume — the KGNO trace has audio to it as well.

Figure 4: Animation of the July 22nd, 2020 earthquake simulation generated by SPECFEM3D Globe. Colour scale represents displacement in centimetres. The left panel shows the unfiltered vertical synthetic seismogram at KGNO, the seismic station located in Kingston, Ontario. Audio tone corresponds to the ups and downs of this trace.

Conclusion

This qualitative comparison reveals that SF3D is capable of simulating an in-phase seismic signal with a similar relative variation in amplitude (displacement) over time. While the displacements of the synthetic signals are not the same as the measured signals, the difference is likely a result of the seismic velocity model used in the simulation.

In addition to the specific results at KGNO, the study shows the capability of SF3D for assessing seismic risk. If we quantify the displacement accuracy of the simulation across all Canadian stations for a given velocity model, we can start thinking about how to create a seismic risk map based on a number of possible earthquake scenarios. Instead of relying on data from earthquakes that have already happened, geoscientists can use simulations to test hypotheses about future scenarios such as a big earthquake occurring along the San Andreas Fault in California, for example.

This is powerful because it represents a huge potential gain for decision-makers at a relatively low cost, since the earthquakes don’t need to happen first. As avid players of SimCity know well, preparation doesn’t count if it only happens after the damage and destruction of the earthquake.

Acknowledgements

Thank you to Dr. Hom Nath Gharti for showing me how to use SPECFEM3D Globe and for patiently assisting me with my supercomputer problems. I also want to thank Dr. Alexander Braun and Dr. Georgia Fotopoulos for providing feedback on the animation. Lastly, I would like to thank MITACS for funding this project.

References

Beyreuther, M., Barsch, R., Krischer, L., Megies, T., Behr, Y., & Wassermann, J. (2010). ObsPy: A Python Toolbox for Seismology. Seismological Research Letters, 81(3), 530–533. https://doi.org/10.1785/gssrl.81.3.530

Geological Survey Of Canada. (1989). Canadian National Seismograph Network. https://doi.org/10.7914/SN/CN

Komatitsch, D., Vilotte, J.-P., Tromp, J., Afanasiev, M., Bozdag, E., Charles, J., Chen, M., Goddeke, D., Hjorleifsdottir, V., Labarta, J., Le Goff, N., Le Loher, P., Liu, Q., Maggi, A., Martin, R., McRitchie, D., Messmer, P., Michea, D., Nissen-Meyer, T., … Zhu, H. (2015). SPECFEM3D GLOBE v7.0.0 [software]. Computational Infrastructure for Geodynamics. https://doi.org/NoDOI

United States Geological Survey. (2020, October 3). M 7.8–99 km SSE of Perryville, Alaska. Earthquake Hazards Program. https://earthquake.usgs.gov/earthquakes/eventpage/us7000asvb/executive

--

--

Benjamin Saadia

Geoscientist interested in mineral exploration and engineering problems and the tools we use to solve them.