Geol 455/655 - Lab 3, Grid Dispersion

For the third exercise, we continue exploring the properties of the numerical simulation. The second one is ``grid dispersion,'' and avoiding it affects the dh parameter and the maximum frequency of our simulation.

The Double Finite Difference

finite-difference equations The first equation at right is the definition of the gradient of a 3-dimensional scalar field. The gradient is a vector that points, well, down the steepest gradient. In acoustic seismology we often make use of the acoustic pressure field P(x,y,z,t) that describes sound waves in time and space. Like pressure P, the temperature T(x,y,z,t) is a scalar field, but one we use to describe heat flow. E3D computes elastic waves as a vector displacement field, which is more complicated but uses the same kind of finite-difference solutions we show here for scalar fields.

The second equation defines the Laplacian, which is just a second spatial derivative of the scalar field P(x,y,z) (partial derivatives over the spatial axes only, ignoring the time t axis). The Laplacian appears in our differential equation for heat flow, and it also appears in differential equations for wave propagation, called wave equations.

Now, one way to solve such differential equations is to find a closed-form analytic expression that solves a simplified case of the differential equation. We looked at an analytic solution for the heat flow equation in one dimension. The way to solve these equations without simplifying them is to finite difference them. This works well for fields that are discretized and exist as discrete samples in computer memory. We start with the definition of the analytic derivative, the third equation at right.

The fourth equation simply removes the limit and keeps dx finite, or greater than zero. The bold Dx is an operator, one that approximates the analytic partial-derivative operator. It says ``apply the single finite difference in x to the field P(x)'' that may vary in the x direction.

The fifth equation shows how to derive the second finite-difference operator Dx2, which is just applying Dx twice to take the difference of the difference and approximate the second spatial derivative.

The sixth line is the result, defining the second-difference operator Dx2 that we can just plug into differential equations containing the Laplacian, along with the first difference, to solve them explicitly. Finite-difference computer programs just loop along time t, in steps dt, to build up a computed wavefield (or temperature field) for the time range we desire. As you saw in lab 2, if we make dt too large in an attempt to save computation time, the faster velocities make the discretized wavefield aliased in time, and the program blows up.

Accuracy of the Double Difference

Double-difference error demo The Courant condition gives us a maximum dt from the maximum velocity in the simulation. The error we can run into using the double difference operator Dx2 to estimate the second spatial derivative will give us a constraint on E3D's dh parameter. This constraint will also depend on the velocities being used in the model being computed. (In E3D, dh = dx = dy = dz; the spatial discretization is the same in all 3 directions.)

At right is a blue sinusoid of amplitude one and wavelength lambda equal to 4 in x. On the left it is sampled with the green line at dx=1 for 4 samples per cycle (or wavelength). On the right it is sampled with the red line at dx=2 for 2 samples per cycle. The coarser sampling and the red line are not aliased, since the Nyquist Criterion says 2 samples per cycle perfectly represent the wave.

The equations on the lower right of the figure here show that the second derivative of this sinusoid is a scaled sinusoid itself. At a value of x equal to three-quarters of the wavelength lambda, which is x=3, the analytic 2nd derivative has a value of 2.47.

Now, let's apply the formula for the 2nd-difference approximation at the same point, x=3. The estimate of the 2nd derivative withDx2 = [P(x+dx)-2P(x)+P(x-dx)]/dx2 is centered at a value of x, which will be 3.

On the left, the 3rd sample on the green line is at a minimum in the sine function. The 2nd derivative is supposed to be about 2.47, and Dx2 produces an estimate of 2.0. For geophysical work, this is sometimes accurate enough.

The second point on the red line, sampling the sine at 2 points per cycle, is at an identical minimum in the sine wave, and ought to still give a 2nd derivative close to 2.47. However, this coarser sampling does not produce nearly as accurate a Dx2 estimate, coming back with a value of 1. The 2nd difference estimate is highly inaccurate, despite the sampling being sufficient for the Nyquist Criterion.

Note that computing Dx2 at x=2, at the sine wave's zero crossing, will always yield the correct value, zero, for the 2nd derivative there. This is true for both the green and red samplings. In fact, Dx2 is near-perfect where the 2nd derivative is small, and underestimates the 2nd derivative more and more at it gets larger.

This factor of two error at the maxima will have drastic effects on the wave propagation in E3D. Waves with high frequencies, having the shortest wavelengths, will seem to slow down and even stop propagating once they get below six spatial samples per wavelength. This slowing of high-frequency parts of the wavefield is reminiscent of the dispersion (spread) of the velocities of surface waves according to their frequency. The computational artifact is thus known as grid dispersion, since it depends on the spatial sampling dh of the model grid.

visual with no grid dispersion visual with grid dispersion For reliable computations, Shawn Larsen and other authors of finite-difference wave modeling codes suggest that the shortest wavelengths in the calculation cover at least 10 grid points-- in other words, that 10dh < lambdamin . At left is an example E3D computation visual (a map with m=0) that has at least 10 spatial samples per wavelength. At right is a visual from the same model, where the source was set to pump a higher frequency into the grid that results in fewer than 6 points per cycle. The grid dispersion generates false long durations of shaking. You can also notice a ``moire'' pattern of waves that seem to be propagating parallel to the axes rather than radiating circularly from the source. As the computation continues, you will see waves propagating so slowly they appear to be stationary.


  1. Given the criterion for avoiding grid dispersion, 10dh < lambdamin , derive an expression that gives you the dh needed to meet this criterion for a highest desired frequency fmax. This expression will involve a velocity-- is this the minimum or the maximum velocity found in the model?

  2. If one had a code that was computing only P-wave propagation in a crustal model, with a dh that avoided grid dispersion, and switched to computing only shear-wave propagation across the same realistic model, would dh have to increase or decrease? Could you guess by how much?

  3. Briefly describe a few examples of lithospheric geologic features that affect the Courant Condition on dt, and then describe the types of features that affect the need to set dh small enough to avoid grid dispersion.

  4. To compute a seismogram representing a given duration, the cost of computation goes up as dt gets smaller, since more time steps are needed. Each time step in general requires computation effort proportional to the number of grid nodes in the model, which is nx times ny times nz. Thus the computational effort required for a model of a certain size goes up as dh decreases, as well. Can you think of a set or sets of lithospheric features, which, if included in a model, cause the computational cost to skyrocket by demanding both small dt and small dh?

  5. Download the E3D input file. Examine it in a text editor. Note that we are not using ModelAssembler for this lab, so all velocities are set in . Given the frequency specified on the source line of the input file, and the other settings in the file, what is the minimum number of grid samples per wavelength for the computation set up here?

    Run this model with the command:

    	/quake/rg/e3d/bin/e3d &
    Watch the visual to verify that no obvious grid dispersion occurs. As you did in Lab 2, rename the LSM and UNLV x-component output files to save them.

  6. Edit the source line of the file to increase the source frequency by a factor of five. How many points per wavelength will this computation have as a minimum? Run this model and observe the grid dispersion artifacts, which become more obvious as the computation proceeds. Again rename the LSM and UNLV output files to save them.

  7. Plot the four seismogram files for the two runs on the LSM and UNLV stations in SAC, label them completely, including which are from which station and which run, and point out some possible features of grid dispersion on the plot.

  8. The original source line in is not a efficient as it might be. How high can you raise the frequency and still avoid grid dispersion? Set this value and make another run. Print out a comparison of the seismograms from these two runs (both without grid dispersion).