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 (psi*psi + eta*eta + q*q);
      double X = sqrt (psi*psi + q*q);
      double dtilde = eta*sind - q*cosd;
      double ytilde = eta*cosd + q*sind;
      double atanp = fabs (q) > 1e-6 ? atan (psi*eta/(q*R)) : 0.;
    
      mulambda = mulambda/(1. + mulambda);
      double logReta = R + eta > 1e-6 ? log (R + eta) : - log (R - eta);
      double Reta = fabs (R + eta) > 1e-6 ? R + eta : 1e30;
      double I1, I2, I3, I4, I5;
      if (fabs (cosd) > 1e-6) {
        /* formula (28) */
        I5 = fabs (psi) < 1e-6 ? 0. :
          mulambda*2./cosd*atan ((eta*(X + q*cosd) + 
    			      X*(R + X)*sind)/(psi*(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*psi/(R + dtilde)) - sind/cosd*I5;
      }
      else {
        /* formula (29) */
        double R1 = R + dtilde;
        I1 = - mulambda/2.*psi*q/(R1*R1);
        I3 = mulambda/2.*(eta/R1 + ytilde*q/(R1*R1) - logReta);
        I2 = mulambda*(- logReta) - I3;
        I4 = - mulambda*q/R1;
        I5 = - mulambda*psi*sind/R1;
      }
        
      /* strike-slip, formula (25) */  
      if (U[0] != 0.) {
        double U1pi = U[0]/(2.*M_PI);
        u[0] -= U1pi*(psi*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 + psi)) + cosd*atanp - I1*sind*cosd);
        u[2] -= U2pi*(dtilde*q/(R*(R + psi)) + 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 + psi)) - 
    		  sind*(psi*q/(R*Reta) - atanp) - I1*sind*sind);
        u[2] += U3pi*(ytilde*q/(R*(R + psi)) + 
    		  cosd*(psi*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 (delta), sind = sin (delta);
      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;
    }
    
    typedef struct {
      double depth, x, y;
      double strike, dip, rake;
      double length, width, U;
      double vU[3];
    } Fault;
    
    struct Okada {
      scalar d;
      double x, y, depth;
      double strike, dip, rake;
      double mu, lambda;
      double length, width, vU[3], U;
      double R;
      int (* iterate) (void);
      bool flat, centroid;
      Fault * faults;
    };
    
    void okada (struct Okada p)
    {
      // default settings
      if (p.mu == 0.)     p.mu = 1.;
      if (p.lambda == 0.) p.lambda = 1.;
      if (p.R == 0.)      p.R = 6371220.; /* Earth radius (metres) */
    
      Fault faults[2] = {0};
      if (p.faults == NULL) {
        faults[0] = (Fault) {p.depth, p.x, p.y,
    			 p.strike, p.dip, p.rake,
    			 p.length, p.width, p.U};
        p.faults = faults;
      }
      foreach() {
        val(p.d) = 0.;
        for (Fault * f = p.faults; f && f->depth > 0.; f++) {
          double depth = f->depth, dtr = pi/180.;
          if (p.centroid)
    	depth -= f->width*fabs (sin (f->dip*dtr))/2.;
          if (f->rake != nodata) {
    	f->vU[0] = f->U*cos (f->rake*dtr);
    	f->vU[1] = f->U*sin (f->rake*dtr);
          }
          double sina = sin ((90. - f->strike)*dtr);
          double cosa = cos ((90. - f->strike)*dtr);
          double sind = sin (f->dip*dtr);
          /* depth of Okada origin */
          depth = sind > 0. ? depth + f->width*sind : depth;
          /* origin to the centroid */
          double x0 = f->length/2., y0 = f->width/2.*cos (f->dip*dtr);
    
          double xc, yc;
          if (p.flat) {
    	xc = x - f->x;
    	yc = y - f->y;
          }
          else {
    	xc = p.R*cos(y*dtr)*dtheta(x, f->x)*dtr;
    	yc = p.R*dtheta(y, f->y)*dtr;
          }
          double x1 =   cosa*xc + sina*yc;
          double y1 = - sina*xc + cosa*yc;
          double oka[3];
          okada_rectangular_source (f->vU, f->length, f->width, depth, 
    				f->dip*dtr,
    				p.mu/p.lambda,
    				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.
    • faults: if non-NULL, this defines an array of several fault parameters which will be used instead of the parameters for a single fault above. This is useful to efficiently define a deformation composed of many Okada subfaults. Note that the array must be terminated by a “dummy fault” of depth smaller than or equal to zero.
    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[];
          eta[] = 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