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;
    }
    
    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;
    };
    
    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) */
    
      double dtr = pi/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.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.
    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