Estimate gravity changes caused by an earthquake

This site interactivelly calculates co- and post-seismic gravity changes of a spherical symmetric, self-gravitating viscoelastic Earth, using a viscosity profile and fault parameters given by a user. The estimated results can be used to interpret satellite gravity field variations. For the detailes, refer to Tanaka et al., 2014. Spectral-finite element approach to post-seismic relaxation in a spherical compressible Earth: application to the gravity field variations due to the 2004 Sumatra-Andaman earthquake, Geophys. J. Int., in press.

How to use this application program?
2 (or 3) text files and a few parameters:
A) a viscosity structure file, B) a source-parameter file, and maximum time step and a cut-off spherical degree (<=100)
(option) C) a site location (longitude & latitude) file, and the radius of a Gaussian filter

D) Spherical harmonic coefficients for the vertical displacement and the potential change,
(if the input C) is given) E) Vertical displacement, geoid height and gravity changes at a given site location

Examples of the input files
5701.E+03 6071.E+03 1.E+19 1.E+19
6071.E+03 6321.E+03 1.E+19 1.E+19
6321.E+03 6371.E+03 1.E+24 1.E+24

Radius 1 (km), Radius 2 (km), Viscosity for Maxwell element (Pas) from radius 1 to radius 2, Viscosity for Kelvin-Voigt element (Pas) from radius 1 to radius 2
Radius 1 < Radius 2, and Radius 1 on the first line and Radius 2 on the last line must be equal to 5701 and 6371 (km), respectively. Below 5701 km, both viscosity values are set at 1.E+21.
In the above example, the lithospheric thickness is 50 km. The format is arbitrary. The code may not work properly for too small viscosity values (e.g. < 1.E+17) because the time step is fixed at 0.1 yr.
0.63680E+07 0.22100E+01 0.94790E+02 0.32200E+03 0.11000E+02 0.91050E+02 0.10000E+01 0.47400E+09
0.63680E+07 0.24400E+01 0.94620E+02 0.32200E+03 0.11000E+02 0.92300E+02 0.10000E+01 0.59200E+09
0.63680E+07 0.26600E+01 0.94440E+02 0.32200E+03 0.11000E+02 0.90440E+02 0.10000E+01 0.54700E+09

1st line: 1st point source, 2nd line: 2nd point source, ...
Each line contains radial distance up to a source [km], lat & lon [deg], strike, dip & rake [deg], slip [m], area [m^2] with an arbitrary format. The radius of the Earth is set to 6371 km in all calculations. The uppermost layer for depths shallower than 3 km is the ocean (PREM). So, place sources at radius <= 6368 km. The maximum number of sources is 1,000 and source depths from 3 km to 100 km are currently available.

C) (option)
80.5400 -10.2300
81.0000 -11.5230
86.2300 -12.1234
The file must be a text file with a file size of up to 2 MB (about 40,000 sites). 1st line: lon & lat [deg] of 1st observation site, 2nd line: lon & lat [deg] of 2nd observation site, ... with an arbitrary format.

Examples of the output files
The effects of the solid earth deformation and the ocean mass redistribution effects can be separately downloaded.

1 0 0.00000000E+00 0.21916190E-03 0.00000000E+00 0.12455367E-02 0.00000000E+00
1 1 0.00000000E+00 0.85570982E-04 0.11538086E-02 -0.86108445E-04 0.38754825E-02
2 0 0.00000000E+00 -0.50411306E-03 0.00000000E+00 -0.44354502E-02 0.00000000E+00
2 1 0.00000000E+00 -0.99001218E-04 0.19046869E-03 0.11833710E-02 0.13643883E-02
2 2 0.00000000E+00 -0.79731904E-03 0.85005452E-04 -0.39899383E-02 0.44096136E-03
3 0 0.00000000E+00 -0.42835218E-03 0.00000000E+00 -0.83018827E-03 0.00000000E+00
3 1 0.00000000E+00 -0.26383104E-04 -0.85190580E-03 -0.60724297E-04 -0.15376649E-02

Each line contains degree j, order m, time [yr], real and imaginary parts of the vertical displacement, Uj,m, and real and imaginary parts of the gravity potential change, Fj,m.

The definition of the spherical harmonic coefficients is described in Martinec, 2000. Spectral-finite element approach to three-dimensional viscoelastic relaxation in a spherical earth, Geophys. J. Int., 142, 117-141.

0.0 80.540 -10.230 -0.31385E-02 -0.91938E-03 -0.52523E+00
0.0 81.000 -11.523 -0.40880E-02 -0.89904E-03 -0.51670E+00
0.0 86.230 -12.123 -0.25478E-02 -0.83130E-03 -0.40576E+00
0.1 80.540 -10.230 -0.33094E-02 -0.93469E-03 -0.55817E+00
0.1 81.000 -11.523 -0.42451E-02 -0.91461E-03 -0.54928E+00
0.1 86.230 -12.123 -0.26433E-02 -0.83181E-03 -0.42066E+00
0.2 80.540 -10.230 -0.34798E-02 -0.95037E-03 -0.59037E+00
0.2 81.000 -11.523 -0.44029E-02 -0.93053E-03 -0.58110E+00
Each line contains time [yr], lon & lat [deg], u_ud [m], geoid height change [m], gravity change [microGal].
The vertical displacement is low-pass filtered, which cannot be directly compared with GNSS displacement data. The geoid height change and the gravity change can be compared with satellite gravity data. It is recommended to use the same cut-off degree and the same radius of the Gaussian filter as used in the analysis of the observed data.

Proceed to computation

Note: This program is a test version and does not gurantee results for any input parameters. When using the results, please cite the above paper (Tanaka et al., 2014).

Update information
Oct. xx, 2014 A test version released. As the density and elastic structure, a 5-layer PREM in the paper but with the top and bottom layers being replaced by the ocean of thickness 3 km and incompressible liquid core is assumed. To save computation time, the number of elements for deeper portions are reduced, which may cause lower accuracy for estimating lower-degree deformations.

Contact: Yoshiyuki Tanaka
Earthquake Research Institute, University of Tokyo
y-tanaka at

Copyright (C) 2014 YT, HT, VK, ZM All Rights Reserved.