src/okada.h

Okada fault model

This is an implementation of the formulae of Okada, 1985.

/* formulae (25)-(30) */
static void rectangular_source (const double U[3], double cosd, double sind,
				double mulambda, double d,
				double psi, double eta, double q,
				double u[3])
{
  double R = sqrt (ψ*ψ + η*η + q*q);
  double X = sqrt (ψ*ψ + q*q);
  double dtilde = η*sind - q*cosd;
  double ytilde = η*cosd + q*sind;
  double atanp = fabs (q) > 1e-6 ? atan (ψ*η/(q*R)) : 0.;

  mulambda = mulambda/(1. + mulambda);
  double logReta = R + η > 1e-6 ? log (R + η) : - log (R - η);
  double Reta = fabs (R + η) > 1e-6 ? R + η : 1e30;
  double I1, I2, I3, I4, I5;
  if (fabs (cosd) > 1e-6) {
    /* formula (28) */
    I5 = fabs (ψ) < 1e-6 ? 0. :
      mulambda*2./cosd*atan ((η*(X + q*cosd) + 
			      X*(R + X)*sind)/(ψ*(R + X)*cosd));
    I4 = mulambda/cosd*(log (R + dtilde) - sind*logReta);
    I3 = mulambda*(1./cosd*ytilde/(R + dtilde) - logReta) + sind/cosd*I4;
    I2 = mulambda*(- logReta) - I3;
    I1 = mulambda*(-1./cosd*ψ/(R + dtilde)) - sind/cosd*I5;
  }
  else {
    /* formula (29) */
    double R1 = R + dtilde;
    I1 = - mulambda/2.*ψ*q/(R1*R1);
    I3 = mulambda/2.*(η/R1 + ytilde*q/(R1*R1) - logReta);
    I2 = mulambda*(- logReta) - I3;
    I4 = - mulambda*q/R1;
    I5 = - mulambda*ψ*sind/R1;
  }
    
  /* strike-slip, formula (25) */  
  if (U[0] != 0.) {
    double U1pi = U[0]/(2.*M_PI);
    u[0] -= U1pi*(ψ*q/(R*Reta) + atanp + I1*sind);
    u[1] -= U1pi*(ytilde*q/(R*Reta) + q*cosd/Reta + I2*sind);
    u[2] -= U1pi*(dtilde*q/(R*Reta) + q*sind/Reta + I4*sind);
  }

  /* dip-slip, formula (26) */  
  if (U[1] != 0.) {
    double U2pi = U[1]/(2.*M_PI);
    u[0] -= U2pi*(q/R - I3*sind*cosd);
    u[1] -= U2pi*(ytilde*q/(R*(R + ψ)) + cosd*atanp - I1*sind*cosd);
    u[2] -= U2pi*(dtilde*q/(R*(R + ψ)) + sind*atanp - I5*sind*cosd);
  }

  /* tensile, formula (27) */  
  if (U[2] != 0.) {
    double U3pi = U[2]/(2.*M_PI);
    u[0] += U3pi*(q*q/(R*Reta) - I3*sind*sind);
    u[1] += U3pi*(-dtilde*q/(R*(R + ψ)) - 
		  sind*(ψ*q/(R*Reta) - atanp) - I1*sind*sind);
    u[2] += U3pi*(ytilde*q/(R*(R + ψ)) + 
		  cosd*(ψ*q/(R*Reta) - atanp) - I5*sind*sind);
  }
}

/* formula (24) */
static void okada_rectangular_source (const double U[3], 
				      double L, double W, double d, 
				      double delta, double mulambda,
				      double x, double y,
				      double u[3])
{
  double cosd = cos (δ), sind = sin (δ);
  double p = y*cosd + d*sind;
  double q = y*sind - d*cosd;

  u[0] = u[1] = u[2] = 0.;
  rectangular_source (U, cosd, sind, mulambda, d,
		      x, p, q,
		      u);
  rectangular_source (U, cosd, sind, mulambda, d,
		      x - L, p - W, q,
		      u);

  double u1[3] = {0., 0., 0.};
  rectangular_source (U, cosd, sind, mulambda, d,
		      x, p - W, q,
		      u1);
  rectangular_source (U, cosd, sind, mulambda, d,
		      x - L, p, q,
		      u1);
  u[0] -= u1[0];
  u[1] -= u1[1];
  u[2] -= u1[2];
}

static double dtheta (double theta1, double theta2)
{
  double d = theta1 - theta2;
  if (d > 180.) d -= 360.;
  if (d < -180.) d += 360.;
  return d;
}

struct Okada {
  scalar d;
  double x, y, depth;
  double strike, dip, rake;
  double μ, λ;
  double length, width, vU[3], U;
  double R;
  int (* iterate) (void);
  bool flat, centroid;
};

void okada (struct Okada p)
{
  // default settings
  if (p.μ == 0.)     p.μ = 1.;
  if (p.λ == 0.) p.λ = 1.;
  if (p.R == 0.)      p.R = 6371220.; /* Earth radius (metres) */

  double dtr = π/180.;
  if (p.centroid)
    p.depth -= p.width*fabs (sin (p.dip*dtr))/2.;
  if (p.rake != nodata) {
    p.vU[0] = p.U*cos (p.rake*dtr);
    p.vU[1] = p.U*sin (p.rake*dtr);
  }
  double sina = sin ((90. - p.strike)*dtr);
  double cosa = cos ((90. - p.strike)*dtr);
  double sind = sin (p.dip*dtr);
  /* depth of Okada origin */
  double depth = sind > 0. ? p.depth + p.width*sind : p.depth;
  /* origin to the centroid */
  double x0 = p.length/2., y0 = p.width/2.*cos (p.dip*dtr);
  
  foreach() {
    if (p.flat) {
      x -= p.x;
      y -= p.y;
    }
    else {
      x = p.R*cos(y*dtr)*dtheta(x, p.x)*dtr;
      y = p.R*dtheta(y, p.y)*dtr;
    }
    double x1 =   cosa*x + sina*y;
    double y1 = - sina*x + cosa*y;
    double oka[3];
    okada_rectangular_source (p.vU, p.length, p.width, depth, 
			      p.dip*dtr,
			      p.μ/p.λ,
			      x0 + x1, y0 + y1,
			      oka);
    val(p.d) = oka[2];
  }
}

User interface

Use function fault() to alter water depth h (where h > dry) according to the fault parameters:

  • x, y: coordinates of the fault centroid (see boolean flat for coordinate type).
  • depth: depth of the top edge of the fault (see also centroid).
  • strike, dip, rake: fault parameters in degrees. (0 <= strike < 360, -90 <= dip <= 90, -90 <= rake <= 90 where rake = 90 degs and dip > 0 is reverse faulting. NB: Okada defines normal faulting by rake = 90 deg and dip < 0 whereas the seismological convention now is generally 0 <= dip <= 90 and rake = -90 for normal faulting).
  • mu, lambda: only the ratio is important and default is mu/lambda = 1.
  • length, width, U: length and width of the fault plane and slip on the fault plane (generally in meters).
  • vU[3]: displacement vector. Note that vU[0] and vU[1] can be calculated from rake and slip (U), vU[2] is tensile opening of a fault so default is 0. (not needed as long as U and rake are given).
  • R: is the radius of the earth (for when x, y are in longitude and latitude i.e. flat = false).
  • iterate: is the function to use to iterate.
  • flat: true assumes x, y cartesian, false assumes longitude and latitude (default).
  • centroid: assumes that depth is measured to the centroid of the fault, not the top edge.
void fault (struct Okada p)
{
  scalar hold[];
  // save the initial water depth
  scalar_clone (hold, h);
  foreach()
    hold[] = h[];
  boundary ({hold});

  p.d = h;
  int nitermax = 20;
  do {
    okada (p);
    // h[] now contains the Okada vertical displacement
    foreach() {
      // deformation is added to hold[] (water depth) only in wet areas
      h[] = hold[] > dry ? max (0., hold[] + h[]) : hold[];
      η[] = zb[] + h[];
    }
  } while (p.iterate && p.iterate() && nitermax--);
}

References

[okada1985]

Yoshimitsu Okada. Surface deformation due to shear and tensile faults in a half-space. Bulletin of the seismological society of America, 75(4):1135-1154, 1985.

Usage

Examples