Smith and Sandwell 
We wish to calculate the displacement vector u(x,y,z) on the surface of the Earth due to a vector body force at depth. This approach is used to describe motion on both curved and discontinuous faults, and is also used to evaluate stress regimes above the upper locking depth. For simplicity, we ignore the effects of EarthÕs sphericity. We assume a Poisson material, and maintain constant moduli with depth. A major difference between this solution and the Okada [1985, 1992] solutions is that we consider deformation due to a vector body force, while the Okada solution considers deformation due to a dislocation. While the following text provides a brief outline of our model formulation, the full derivation and source code of our semi-analytic Fourier model can be found at http://topex.ucsd.edu/body_force/elastic. Our solution is obtained as follows:
(2) Take the 3-D Fourier transform of equations A1 and A2 to reduce the partial differential equations to a set of linear algebraic equations.
(4) Perform the inverse Fourier transform in the z-direction (depth) by repeated application of the Cauchy Residue Theorem. We assume z to be positive upward.
(5) Solve the Boussinesq problem to correct for non-zero normal traction on the half-space. This derivation follows the approach of Steketee  where we impose a negative surface traction in an elastic half space in order to cancel the non-zero traction from the source and image in the elastic full-space.
(6) Integrate the point source Green's function to simulate a fault. For a complex dipping fault, this integration could be done numerically. However, if the faults are assumed to be vertical, the integration can be performed analytically. The body force is applied between the lower depth d1 (e.g., minus infinity) and the upper depth d2 (Figure A1). The displacement or stress (derivatives are computed analytically) can be evaluated at any depth z above d2. Note that the full displacement solution is the sum of three terms: a source, an image, and a Boussinesq correction.
The individual elements of the source and image tensors are
The individual elements of the Boussinesq correction are
(7) Construct a force couple by taking the derivative of the point source in the direction normal to the fault trace. In practice, the body forces due to the stress discontinuity across a fault plane are approximated by the derivative of a Gaussian function, effectively producing a model fault with a finite thickness (Figure A1). Curved faults are constructed with overlapping line segments having cosine tapered ends and are typically 6-10 km long.
The fault trace is imbedded in a two-dimensional grid which is Fourier transformed, multiplied by the transfer functions above (5-7), and inverse Fourier transformed. A constant shear modulus (4.12 x 1010 Pa) and Poisson ratio (0.25) were adopted for all calculations. When the lower edge of the fault is extended to infinite depth, as in the case of the SAF system model, a Fourier cosine transform (mirrored pair) is used in the across-fault direction to maintain the far-field velocity V step across the plate boundary, effectively conserving moment within the grid. Note that this requires the velocity-difference (i.e., stress drop) across a system of connecting faults to have a constant value . To avoid Fourier artifacts where the fault enters the bottom of the grid and leaves the top of the grid, the fault is extended beyond the top of the model and angled to match the intersection point at the bottom. (Figure 1) In addition to computing the velocity field, strain and stress rates are computed from the derivatives of the model. Horizontal derivatives are computed by multiplication of in the Fourier transform domain and vertical derivatives are computed analytically from the transfer functions in A5 and A6.
Equations (4-7) were checked using the computer algebra capabilities in Matlab and then compared to the simple arctangent function for a two-dimensional fault [Weertman, 1964]
Of course, in the two-dimensional case, the Boussinesq correction is not needed and equation (A-5) reduces to the above arctangent formula. Figure 2 provides a numerical comparison between the above arctangent function (analytic profile) and our semi-analytic Fourier solution (Fourier profile). The numerical solutions have relative errors less than 1% as long as the observation depth z is more than one grid cell size above the locking depth d2.
The numerical approach is very efficient; for example, writing and displaying of the deformation/stress grids requires more computer time than the actual computation of the model. The horizontal complexity of the fault system has no effect on the speed of the computation. However, variations in locking depth along the fault system require computing the model for each different locking depth and summing the outputs to form the full solution. The extension to a viscoelastic half space would not introduce a computational burden on an ordinary workstation.