diff --git a/include/materials/linear_elastic.h b/include/materials/linear_elastic.h
index 2ee3271fd..fb5a7eaaf 100644
--- a/include/materials/linear_elastic.h
+++ b/include/materials/linear_elastic.h
@@ -77,6 +77,10 @@ class LinearElastic : public Material<Tdim> {
   double poisson_ratio_{std::numeric_limits<double>::max()};
   //! Bulk modulus
   double bulk_modulus_{std::numeric_limits<double>::max()};
+  //! Compressional Wave Velocity
+  double vp_{std::numeric_limits<double>::max()};
+  //! Shear Wave Velocity
+  double vs_{std::numeric_limits<double>::max()};
 };  // LinearElastic class
 }  // namespace mpm
 
diff --git a/include/materials/linear_elastic.tcc b/include/materials/linear_elastic.tcc
index e1e3d3859..665a1f677 100644
--- a/include/materials/linear_elastic.tcc
+++ b/include/materials/linear_elastic.tcc
@@ -9,10 +9,24 @@ mpm::LinearElastic<Tdim>::LinearElastic(unsigned id,
         material_properties.at("youngs_modulus").template get<double>();
     poisson_ratio_ =
         material_properties.at("poisson_ratio").template get<double>();
+
     // Calculate bulk modulus
     bulk_modulus_ = youngs_modulus_ / (3.0 * (1. - 2. * poisson_ratio_));
 
+    // Calculate constrained and shear modulus
+    double constrained_modulus =
+        youngs_modulus_ * (1. - poisson_ratio_) /
+        ((1. + poisson_ratio_) * (1. - 2. * poisson_ratio_));
+    double shear_modulus = youngs_modulus_ / (2.0 * (1. + poisson_ratio_));
+
+    // Calculate wave velocities
+    vp_ = sqrt(constrained_modulus / density_);
+    vs_ = sqrt(shear_modulus / density_);
+
     properties_ = material_properties;
+    properties_["pwave_velocity"] = vp_;
+    properties_["swave_velocity"] = vs_;
+
     // Set elastic tensor
     this->compute_elastic_tensor();
   } catch (Json::exception& except) {
diff --git a/tests/materials/linear_elastic_test.cc b/tests/materials/linear_elastic_test.cc
index 791cc7de5..75dc611e9 100644
--- a/tests/materials/linear_elastic_test.cc
+++ b/tests/materials/linear_elastic_test.cc
@@ -71,8 +71,6 @@ TEST_CASE("LinearElastic is checked in 2D", "[material][linear_elastic][2D]") {
             Approx(jmaterial["density"]).epsilon(Tolerance));
     REQUIRE(material->template property<double>("youngs_modulus") ==
             Approx(jmaterial["youngs_modulus"]).epsilon(Tolerance));
-    REQUIRE(material->template property<double>("poisson_ratio") ==
-            Approx(jmaterial["poisson_ratio"]).epsilon(Tolerance));
 
     // Check if state variable is initialised
     SECTION("State variable is initialised") {
@@ -147,6 +145,20 @@ TEST_CASE("LinearElastic is checked in 2D", "[material][linear_elastic][2D]") {
     REQUIRE(stress(4) == Approx(0.00000000000000e+00).epsilon(Tolerance));
     REQUIRE(stress(5) == Approx(0.00000000000000e+00).epsilon(Tolerance));
   }
+
+  SECTION("LinearElastic check properties earthquake") {
+    unsigned id = 0;
+
+    auto material =
+        Factory<mpm::Material<Dim>, unsigned, const Json&>::instance()->create(
+            "LinearElastic2D", std::move(id), jmaterial);
+
+    // Get P-Wave and S-Wave Velocities
+    REQUIRE(material->template property<double>("pwave_velocity") ==
+            Approx(116.023870223).epsilon(Tolerance));
+    REQUIRE(material->template property<double>("swave_velocity") ==
+            Approx(62.0173672946).epsilon(Tolerance));
+  }
 }
 
 //! Check linearelastic class in 3D
@@ -289,4 +301,18 @@ TEST_CASE("LinearElastic is checked in 3D", "[material][linear_elastic][3D]") {
     REQUIRE(stress(4) == Approx(7.69230769230769e+01).epsilon(Tolerance));
     REQUIRE(stress(5) == Approx(1.15384615384615e+02).epsilon(Tolerance));
   }
+
+  SECTION("LinearElastic check properties earthquake") {
+    unsigned id = 0;
+
+    auto material =
+        Factory<mpm::Material<Dim>, unsigned, const Json&>::instance()->create(
+            "LinearElastic3D", std::move(id), jmaterial);
+
+    // Get P-Wave and S-Wave Velocities
+    REQUIRE(material->template property<double>("pwave_velocity") ==
+            Approx(116.023870223).epsilon(Tolerance));
+    REQUIRE(material->template property<double>("swave_velocity") ==
+            Approx(62.0173672946).epsilon(Tolerance));
+  }
 }