Merge pull request #134 from kwschultz/master

Fixing unit conventions for moment magnitude

Author: kwschultz

quakelib/src/QuakeLibIO.cpp

@@ -756,7 +756,7 @@ void quakelib::ModelWorld::create_section(std::vector<unsigned int> &unused_trac
     double moment_magnitude = 4.0+log10(conv.sqm2sqkm(fault_area));
 
     for (unsigned int i=0; i<elem_ids.size(); ++i) {
-        double max_slip = pow(10, (3.0/2.0)*(moment_magnitude+10.7))/(1e7*element(elem_ids[i]).lame_mu()*fault_area);
+        double max_slip = pow(10, (3.0/2.0)*(moment_magnitude+6.0))/(element(elem_ids[i]).lame_mu()*fault_area);
         element(elem_ids[i]).set_max_slip(max_slip);
     }
 }
@@ -872,7 +872,8 @@ void quakelib::ModelWorld::create_faults(const std::string &taper_method) {
 // Schultz: Adapted from Steve Ward's model, used in the EQSim comparison for UCERF2 model
 void quakelib::ModelWorld::compute_stress_drops(const double &stress_drop_factor) {
     std::map<UIndex, ModelElement>::iterator eit;
-    UIndex fault_id;
+    UIndex     fault_id;
+    Conversion conv;
     ModelFault this_fault;
     double fault_area, fault_length, fault_width, char_magnitude, char_slip, R, nu, stress_drop;
 
@@ -886,8 +887,8 @@ void quakelib::ModelWorld::compute_stress_drops(const double &stress_drop_factor
         fault_length = this_fault.length();
         fault_width = fault_area/fault_length;
 
-        char_magnitude = 4.0+log10(fault_area*1e-6) + stress_drop_factor;
-        char_slip = pow(10, (3.0/2.0)*(char_magnitude+10.7))/(1e7*eit->second.lame_mu()*fault_area);
+        char_magnitude = 4.0+log10(conv.sqm2sqkm(fault_area)) + stress_drop_factor;
+        char_slip = pow(10, (3.0/2.0)*(char_magnitude+6.0))/(eit->second.lame_mu()*fault_area);
 
         nu = 0.5*eit->second.lame_lambda()/(eit->second.lame_mu() + eit->second.lame_lambda());
         R  = sqrt(fault_width*fault_width + fault_length*fault_length);
@@ -2254,7 +2255,7 @@ int quakelib::ModelWorld::read_files_eqsim(const std::string &geom_file_name, co
         double moment_magnitude = 4.0+log10(conv.sqm2sqkm(fault_areas[eit->section_id()]));
 
         // Need to document where this scaling law comes from
-        double max_slip = pow(10, (3.0/2.0)*(moment_magnitude+10.7))/(1e7*(eit->lame_mu())*fault_areas[eit->section_id()]);
+        double max_slip = pow(10, (3.0/2.0)*(moment_magnitude+6.0))/((eit->lame_mu())*fault_areas[eit->section_id()]);
 
         // Set the max slip for the current element
         eit->set_max_slip(max_slip);

quakelib/src/QuakeLibIO.h

@@ -898,7 +898,7 @@ namespace quakelib {
                     moment += it->_slip*it->_mu*it->_area;
                 }
 
-                _data._event_magnitude = (2.0/3.0)*log10(1e7*moment) - 10.7;
+                _data._event_magnitude = (2.0/3.0)*log10(moment) - 6.0;
             }
 
             //! Get the total amount a given block slipped during this event
@@ -930,7 +930,7 @@ namespace quakelib {
                     moment += getBlockMu(*it)*getEventSlip(*it)*getEventArea(*it);
                 }
 
-                return (2.0/3.0)*log10(1e7*moment) - 10.7;
+                return (2.0/3.0)*log10(moment) - 6.0;
             }
 
             //! Get the magnitude of the earthquake in this event.

src/core/Simulation.h

@@ -233,7 +233,7 @@ class Simulation : public SimFramework, public VCParams, public VCSimData, publi
 
             // Get current expected slip from the current event area
             char_magnitude = 4.0+log10(current_event_area*1e-6) + stressDropFactor();
-            char_slip = pow(10, (3.0/2.0)*(char_magnitude+10.7))/(1e7*getBlock(gid).lame_mu()*current_event_area);
+            char_slip = pow(10, (3.0/2.0)*(char_magnitude+6.0))/(getBlock(gid).lame_mu()*current_event_area);
 
             dynamicStressDrop = -2*getBlock(gid).lame_mu()*char_slip*( (1-nu)*fault_length/fault_width + fault_width/fault_length )/( (1-nu)*M_PI*R );
 

src/simulation/RunEvent.cpp

@@ -811,7 +811,7 @@ void RunEvent::processAftershock(Simulation *sim) {
 
         // Determine the amount of slip needed to match the aftershock magnitude
         // The contribution of each block to moment is based on its fraction of total area*mu
-        double total_moment = pow(10, (as.mag + 10.7)*(3.0/2.0))/1e7;
+        double total_moment = pow(10, (as.mag + 6.0)*(3.0/2.0));
 
         for (bit=id_set.begin(); bit!=id_set.end(); ++bit) {
             Block &b=sim->getBlock(*bit);