diff --git a/include/data_types.h b/include/data_types.h
index 061d27e81..04ea207b9 100644
--- a/include/data_types.h
+++ b/include/data_types.h
@@ -34,6 +34,13 @@ inline double zero() {
   return 0.;
 }
 
+//! Position type
+//! None: No position is specified
+//! Corner: Nodes at boundary corners
+//! Edge: Nodes along boundary edges
+//! Face: Nodes on boundary faces
+enum class Position { None, Corner, Edge, Face };
+
 }  // namespace mpm
 
 #endif  // MPM_DATA_TYPES_H_
diff --git a/include/loads_bcs/absorbing_constraint.h b/include/loads_bcs/absorbing_constraint.h
new file mode 100644
index 000000000..48503c9a1
--- /dev/null
+++ b/include/loads_bcs/absorbing_constraint.h
@@ -0,0 +1,70 @@
+#ifndef MPM_ABSORBING_CONSTRAINT_H_
+#define MPM_ABSORBING_CONSTRAINT_H_
+
+#include "data_types.h"
+
+namespace mpm {
+
+//! AbsorbingConstraint class to store friction constraint on a set
+//! \brief AbsorbingConstraint class to store a constraint on a set
+//! \details AbsorbingConstraint stores the constraint as a static value
+class AbsorbingConstraint {
+ public:
+  // Constructor
+  //! \param[in] setid  set id
+  //! \param[in] dir Direction of p-wave propagation in model
+  //! \param[in] delta Virtual viscous layer thickness
+  //! \param[in] h_min Characteristic length (cell height)
+  //! \param[in] a Dimensionless dashpot weight factor, p-wave
+  //! \param[in] b Dimensionless dashpot weight factor, s-wave
+  //! \param[in] position Nodal position along boundary
+  AbsorbingConstraint(int setid, unsigned dir, double delta, double h_min,
+                      double a = 1., double b = 1.,
+                      mpm::Position position = mpm::Position::None)
+      : setid_{setid},
+        dir_{dir},
+        delta_{delta},
+        h_min_{h_min},
+        a_{a},
+        b_{b},
+        position_{position} {};
+
+  // Set id
+  int setid() const { return setid_; }
+
+  // Direction
+  unsigned dir() const { return dir_; }
+
+  // Virtual viscous layer thickness
+  double delta() const { return delta_; }
+
+  // Cell height
+  double h_min() const { return h_min_; }
+
+  // Return a
+  double a() const { return a_; }
+
+  // Return b
+  double b() const { return b_; }
+
+  // Return position
+  mpm::Position position() const { return position_; }
+
+ private:
+  // ID
+  int setid_;
+  // Direction
+  unsigned dir_;
+  // Delta
+  double delta_;
+  // Cell height
+  double h_min_;
+  // a
+  double a_;
+  // b
+  double b_;
+  // Node position
+  mpm::Position position_;
+};
+}  // namespace mpm
+#endif  // MPM_ABSORBING_CONSTRAINT_H_
diff --git a/include/loads_bcs/constraints.h b/include/loads_bcs/constraints.h
index 0efdef6f6..42c533863 100644
--- a/include/loads_bcs/constraints.h
+++ b/include/loads_bcs/constraints.h
@@ -3,6 +3,7 @@
 
 #include <memory>
 
+#include "absorbing_constraint.h"
 #include "friction_constraint.h"
 #include "logger.h"
 #include "mesh.h"
@@ -48,6 +49,22 @@ class Constraints {
       const std::vector<std::tuple<mpm::Index, unsigned, int, double>>&
           friction_constraints);
 
+  //! Assign nodal absorbing constraints
+  //! \param[in] setid Node set id
+  //! \param[in] absorbing_constraints Constraint at node, dir, delta, h_min, a,
+  //! b
+  bool assign_nodal_absorbing_constraint(
+      int nset_id,
+      const std::shared_ptr<mpm::AbsorbingConstraint>& absorbing_constraints);
+
+  //! Assign absorbing constraints to nodes
+  //! \param[in] absorbing_constraints Constraint at node, dir, delta, h_min, a,
+  //! b, and position
+  bool assign_nodal_absorbing_constraints(
+      const std::vector<std::tuple<mpm::Index, unsigned, double, double, double,
+                                   double, mpm::Position>>&
+          absorbing_constraints);
+
  private:
   //! Mesh object
   std::shared_ptr<mpm::Mesh<Tdim>> mesh_;
diff --git a/include/loads_bcs/constraints.tcc b/include/loads_bcs/constraints.tcc
index 5813e609f..3c66f5130 100644
--- a/include/loads_bcs/constraints.tcc
+++ b/include/loads_bcs/constraints.tcc
@@ -105,3 +105,83 @@ bool mpm::Constraints<Tdim>::assign_nodal_friction_constraints(
   }
   return status;
 }
+
+//! Apply absorbing constraints to nodes
+template <unsigned Tdim>
+bool mpm::Constraints<Tdim>::assign_nodal_absorbing_constraint(
+    int nset_id,
+    const std::shared_ptr<mpm::AbsorbingConstraint>& absorbing_constraint) {
+  bool status = true;
+  try {
+    // int set_id = absorbing_constraint->setid();
+    int set_id = nset_id;
+    auto nset = mesh_->nodes(set_id);
+    if (nset.size() == 0)
+      throw std::runtime_error(
+          "Node set is empty for application of absorbing constraints");
+    unsigned dir = absorbing_constraint->dir();
+    double delta = absorbing_constraint->delta();
+    double h_min = absorbing_constraint->h_min();
+    double a = absorbing_constraint->a();
+    double b = absorbing_constraint->b();
+    mpm::Position position = absorbing_constraint->position();
+    if (delta >= h_min / (2 * a) and delta >= h_min / (2 * b))
+      for (auto nitr = nset.cbegin(); nitr != nset.cend(); ++nitr) {
+        if (!(*nitr)->apply_absorbing_constraint(dir, delta, h_min, a, b,
+                                                 position))
+          throw std::runtime_error(
+              "Failed to apply absorbing constraint at node");
+      }
+    else
+      throw std::runtime_error("Invalid value for delta");
+  } catch (std::exception& exception) {
+    console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
+    status = false;
+  }
+  return status;
+}
+
+//! Assign absorbing constraints to nodes
+template <unsigned Tdim>
+bool mpm::Constraints<Tdim>::assign_nodal_absorbing_constraints(
+    const std::vector<std::tuple<mpm::Index, unsigned, double, double, double,
+                                 double, mpm::Position>>&
+        absorbing_constraints) {
+  bool status = true;
+  try {
+    for (const auto& absorbing_constraint : absorbing_constraints) {
+      // Node id
+      mpm::Index nid = std::get<0>(absorbing_constraint);
+      // Direction
+      unsigned dir = std::get<1>(absorbing_constraint);
+      // Delta
+      double delta = std::get<2>(absorbing_constraint);
+      // h_min
+      double h_min = std::get<3>(absorbing_constraint);
+      // a
+      double a = std::get<4>(absorbing_constraint);
+      // b
+      double b = std::get<5>(absorbing_constraint);
+      // Position
+      mpm::Position position = std::get<6>(absorbing_constraint);
+      // delta check
+      if (delta >= h_min / (2 * a) and delta >= h_min / (2 * b)) {
+        if (position == mpm::Position::Corner or
+            position == mpm::Position::Edge or
+            position == mpm::Position::Face) {
+          // Apply constraint
+          if (!mesh_->node(nid)->apply_absorbing_constraint(dir, delta, h_min,
+                                                            a, b, position))
+            throw std::runtime_error(
+                "Nodal absorbing constraints assignment failed");
+        } else
+          throw std::runtime_error("Invalid position");
+      } else
+        throw std::runtime_error("Invalid value for delta");
+    }
+  } catch (std::exception& exception) {
+    console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
+    status = false;
+  }
+  return status;
+}
\ No newline at end of file
diff --git a/include/mesh.h b/include/mesh.h
index bb51f161a..10193dbda 100644
--- a/include/mesh.h
+++ b/include/mesh.h
@@ -24,6 +24,7 @@
 #include "json.hpp"
 using Json = nlohmann::json;
 
+#include "absorbing_constraint.h"
 #include "cell.h"
 #include "factory.h"
 #include "friction_constraint.h"
diff --git a/include/mesh.tcc b/include/mesh.tcc
index ba0071551..58fecdb51 100644
--- a/include/mesh.tcc
+++ b/include/mesh.tcc
@@ -2000,6 +2000,10 @@ void mpm::Mesh<Tdim>::create_nodal_properties() {
                                        materials_.size());
     nodal_properties_->create_property("normal_unit_vectors", nrows,
                                        materials_.size());
+    nodal_properties_->create_property("wave_velocities", nrows,
+                                       materials_.size());
+    nodal_properties_->create_property("density", nodes_.size(),
+                                       materials_.size());
 
     // Iterate over all nodes to initialise the property handle in each node
     // and assign its node id as the prop id in the nodal property data pool
diff --git a/include/node.h b/include/node.h
index 31a6b9a03..7734d26c0 100644
--- a/include/node.h
+++ b/include/node.h
@@ -213,6 +213,18 @@ class Node : public NodeBase<Tdim> {
   //! \param[in] dt Time-step
   void apply_friction_constraints(double dt) override;
 
+  //! Apply absorbing constraint
+  //! Directions can take values between 0 and Dim * Nphases
+  //! \param[in] dir Direction of p-wave propagation in model
+  //! \param[in] delta Virtual viscous layer thickness
+  //! \param[in] h_min Characteristic length (cell height)
+  //! \param[in] a Dimensionless dashpot weight factor, p-wave
+  //! \param[in] b Dimensionless dashpot weight factor, s-wave
+  //! \param[in] position Nodal position along boundary
+  bool apply_absorbing_constraint(unsigned dir, double delta, double h_min,
+                                  double a, double b,
+                                  mpm::Position position) override;
+
   //! Assign rotation matrix
   //! \param[in] rotation_matrix Rotation matrix of the node
   void assign_rotation_matrix(
diff --git a/include/node.tcc b/include/node.tcc
index 7d3a9cfcb..ec883f928 100644
--- a/include/node.tcc
+++ b/include/node.tcc
@@ -346,6 +346,94 @@ void mpm::Node<Tdim, Tdof, Tnphases>::apply_velocity_constraints() {
   }
 }
 
+// !Apply absorbing constraints
+template <unsigned Tdim, unsigned Tdof, unsigned Tnphases>
+bool mpm::Node<Tdim, Tdof, Tnphases>::apply_absorbing_constraint(
+    unsigned dir, double delta, double h_min, double a, double b,
+    mpm::Position position) {
+  bool status = true;
+  try {
+
+    if (dir >= Tdim) {
+      throw std::runtime_error("Direction is out of bounds");
+      status = false;
+    }
+    if (material_ids_.size() != 0) {
+      // Get material id
+      auto mat_id = material_ids_.begin();
+
+      // Extract material properties and displacements
+      double pwave_v = this->property_handle_->property(
+          "wave_velocities", prop_id_, *mat_id, Tdim)(0);
+      double swave_v = this->property_handle_->property(
+          "wave_velocities", prop_id_, *mat_id, Tdim)(1);
+      double density =
+          this->property_handle_->property("density", prop_id_, *mat_id)(0);
+      Eigen::Matrix<double, Tdim, 1> material_displacement =
+          this->property_handle_->property("displacements", prop_id_, *mat_id,
+                                           Tdim);
+      material_displacement /= this->mass(*mat_id);
+
+      // Wave velocity Eigen Matrix
+      Eigen::Matrix<double, Tdim, 1> wave_velocity =
+          Eigen::MatrixXd::Constant(Tdim, 1, b * swave_v);
+      wave_velocity(dir, 0) = a * pwave_v;
+      // Spring constant Eigen matrix
+      double k_s = (density * pow(swave_v, 2)) / delta;
+      double k_p = (density * pow(pwave_v, 2)) / delta;
+      Eigen::Matrix<double, Tdim, 1> spring_constant =
+          Eigen::MatrixXd::Constant(Tdim, 1, k_s);
+      spring_constant(dir, 0) = k_p;
+
+      // Iterate through each phase
+      for (unsigned phase = 0; phase < Tnphases; ++phase) {
+        // Calculate Aborbing Traction
+        Eigen::Matrix<double, Tdim, 1> absorbing_traction_ =
+            this->velocity_.col(phase).cwiseProduct(wave_velocity) * density +
+            material_displacement.cwiseProduct(spring_constant);
+
+        // Update external force
+        if (Tdim == 2) {
+          Eigen::Matrix<double, Tdim, 1> absorbing_force_;
+          switch (position) {
+            case mpm::Position::Corner:
+              absorbing_force_ = 0.5 * h_min * absorbing_traction_;
+              break;
+            case mpm::Position::Edge:
+              absorbing_force_ = h_min * absorbing_traction_;
+              break;
+            default:
+              throw std::runtime_error("Invalid absorbing boundary position");
+              break;
+          }
+          this->update_external_force(true, phase, -absorbing_force_);
+        } else if (Tdim == 3) {
+          Eigen::Matrix<double, Tdim, 1> absorbing_force_;
+          switch (position) {
+            case mpm::Position::Corner:
+              absorbing_force_ = 0.25 * pow(h_min, 2) * absorbing_traction_;
+              break;
+            case mpm::Position::Edge:
+              absorbing_force_ = 0.5 * pow(h_min, 2) * absorbing_traction_;
+              break;
+            case mpm::Position::Face:
+              absorbing_force_ = pow(h_min, 2) * absorbing_traction_;
+              break;
+            default:
+              throw std::runtime_error("Invalid absorbing boundary position");
+              break;
+          }
+          this->update_external_force(true, phase, -absorbing_force_);
+        }
+      }
+    }
+  } catch (std::exception& exception) {
+    console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
+    status = false;
+  }
+  return status;
+}
+
 //! Assign friction constraint
 //! Constrain directions can take values between 0 and Dim * Nphases
 template <unsigned Tdim, unsigned Tdof, unsigned Tnphases>
@@ -543,10 +631,15 @@ void mpm::Node<Tdim, Tdof, Tnphases>::update_property(
     bool update, const std::string& property,
     const Eigen::MatrixXd& property_value, unsigned mat_id,
     unsigned nprops) noexcept {
+
   // Update/assign property
   node_mutex_.lock();
-  property_handle_->update_property(property, prop_id_, mat_id, property_value,
-                                    nprops);
+  if (update)
+    property_handle_->update_property(property, prop_id_, mat_id,
+                                      property_value, nprops);
+  else
+    property_handle_->assign_property(property, prop_id_, mat_id,
+                                      property_value, nprops);
   node_mutex_.unlock();
 }
 
diff --git a/include/node_base.h b/include/node_base.h
index 17eddfdca..f585e8bc0 100644
--- a/include/node_base.h
+++ b/include/node_base.h
@@ -205,6 +205,18 @@ class NodeBase {
   //! \param[in] dt Time-step
   virtual void apply_friction_constraints(double dt) = 0;
 
+  //! Apply absorbing constraint
+  //! Directions can take values between 0 and Dim * Nphases
+  //! \param[in] dir Direction of p-wave propagation in model
+  //! \param[in] delta Virtual viscous layer thickness
+  //! \param[in] h_min Characteristic length (cell height)
+  //! \param[in] a Dimensionless dashpot weight factor, p-wave
+  //! \param[in] b Dimensionless dashpot weight factor, s-wave
+  //! \param[in] position Nodal position along boundary
+  virtual bool apply_absorbing_constraint(unsigned dir, double delta,
+                                          double h_min, double a, double b,
+                                          mpm::Position position) = 0;
+
   //! Assign rotation matrix
   //! \param[in] rotation_matrix Rotation matrix of the node
   virtual void assign_rotation_matrix(
diff --git a/include/particles/particle.h b/include/particles/particle.h
index a362d7c76..f1dc40d5d 100644
--- a/include/particles/particle.h
+++ b/include/particles/particle.h
@@ -149,6 +149,9 @@ class Particle : public ParticleBase<Tdim> {
   //! Map multimaterial domain gradients to nodes
   void map_multimaterial_domain_gradients_to_nodes() noexcept override;
 
+  // ! Map linear elastic wave velocities to nodes
+  void map_wave_velocities_to_nodes() noexcept override;
+
   //! Assign nodal mass to particles
   //! \param[in] mass Mass from the particles in a cell
   //! \retval status Assignment status
diff --git a/include/particles/particle.tcc b/include/particles/particle.tcc
index b50a5eb56..0f3203e6d 100644
--- a/include/particles/particle.tcc
+++ b/include/particles/particle.tcc
@@ -582,6 +582,30 @@ void mpm::Particle<
   }
 }
 
+//! Map linear elastic wave velocities to nodes
+template <unsigned Tdim>
+void mpm::Particle<Tdim>::map_wave_velocities_to_nodes() noexcept {
+  // if (std::isnan(this->material())->template
+  // property<double>(std::string("pwave_velocity")) == false) {
+  Eigen::Matrix<double, 1, 1> density;
+  density(0) =
+      (this->material())->template property<double>(std::string("density"));
+  Eigen::Matrix<double, Tdim, 1> wave_velocities;
+  wave_velocities(0) =
+      (this->material())
+          ->template property<double>(std::string("pwave_velocity"));
+  wave_velocities(1) =
+      (this->material())
+          ->template property<double>(std::string("swave_velocity"));
+  for (unsigned i = 0; i < nodes_.size(); ++i) {
+    nodes_[i]->update_property(false, "wave_velocities", wave_velocities,
+                               this->material_id(), Tdim);
+    nodes_[i]->update_property(false, "density", density, this->material_id(),
+                               1);
+    //}
+  }
+}
+
 // Compute strain rate of the particle
 template <>
 inline Eigen::Matrix<double, 6, 1> mpm::Particle<1>::compute_strain_rate(
diff --git a/include/particles/particle_base.h b/include/particles/particle_base.h
index dcb2bc419..260af1eac 100644
--- a/include/particles/particle_base.h
+++ b/include/particles/particle_base.h
@@ -159,6 +159,9 @@ class ParticleBase {
   //! Map multimaterial domain gradients to nodes
   virtual void map_multimaterial_domain_gradients_to_nodes() noexcept = 0;
 
+  // ! Map linear elastic wave velocities to nodes
+  virtual void map_wave_velocities_to_nodes() noexcept = 0;
+
   //! Assign material
   virtual bool assign_material(const std::shared_ptr<Material<Tdim>>& material,
                                unsigned phase = mpm::ParticlePhase::Solid) = 0;
diff --git a/include/solvers/mpm_base.h b/include/solvers/mpm_base.h
index c84eaf1f8..52f6b6575 100644
--- a/include/solvers/mpm_base.h
+++ b/include/solvers/mpm_base.h
@@ -99,6 +99,9 @@ class MPMBase : public MPM {
   //! Particle velocity constraints
   void particle_velocity_constraints();
 
+  //! Apply Absorbing Constraints
+  void nodal_absorbing_constraints();
+
  private:
   //! Return if a mesh will be isoparametric or not
   //! \retval isoparametric Status of mesh type
@@ -127,6 +130,12 @@ class MPMBase : public MPM {
   void nodal_frictional_constraints(
       const Json& mesh_prop, const std::shared_ptr<mpm::IOMesh<Tdim>>& mesh_io);
 
+  //! Nodal absorbing constraints
+  //! \param[in] mesh_prop Mesh properties
+  //! \param[in] mesh_io Mesh IO handle
+  void nodal_absorbing_constraints(
+      const Json& mesh_prop, const std::shared_ptr<mpm::IOMesh<Tdim>>& mesh_io);
+
   //! Cell entity sets
   //! \param[in] mesh_prop Mesh properties
   //! \param[in] check Check duplicates
@@ -207,6 +216,11 @@ class MPMBase : public MPM {
   double damping_factor_{0.};
   //! Locate particles
   bool locate_particles_{true};
+  //! Absorbing Boundary Variables
+  bool absorbing_boundary_{false};
+  std::vector<std::shared_ptr<mpm::AbsorbingConstraint>> absorbing_constraint_;
+  std::vector<int> absorbing_nset_id_;
+  mpm::Position position_{mpm::Position::None};
 
 #ifdef USE_GRAPH_PARTITIONING
   // graph pass the address of the container of cell
diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc
index f829b5f95..d71d53df9 100644
--- a/include/solvers/mpm_base.tcc
+++ b/include/solvers/mpm_base.tcc
@@ -236,6 +236,9 @@ void mpm::MPMBase<Tdim>::initialise_mesh() {
   // Read and assign friction constraints
   this->nodal_frictional_constraints(mesh_props, mesh_io);
 
+  // Read and assign absorbing constraintes
+  this->nodal_absorbing_constraints(mesh_props, mesh_io);
+
   // Initialise cell
   auto cells_begin = std::chrono::steady_clock::now();
   // Shape function name
@@ -937,6 +940,78 @@ void mpm::MPMBase<Tdim>::nodal_frictional_constraints(
   }
 }
 
+// Assign nodal absorbing constraints
+template <unsigned Tdim>
+void mpm::MPMBase<Tdim>::nodal_absorbing_constraints(
+    const Json& mesh_props, const std::shared_ptr<mpm::IOMesh<Tdim>>& mesh_io) {
+  try {
+    // Read and assign absorbing constraints
+    if (mesh_props.find("boundary_conditions") != mesh_props.end() &&
+        mesh_props["boundary_conditions"].find("absorbing_constraints") !=
+            mesh_props["boundary_conditions"].end()) {
+      // Iterate over absorbing constraints
+      for (const auto& constraints :
+           mesh_props["boundary_conditions"]["absorbing_constraints"]) {
+        // Set id
+        int nset_id = constraints.at("nset_id").template get<int>();
+        // Direction
+        unsigned dir = constraints.at("dir").template get<unsigned>();
+        // Delta
+        double delta = constraints.at("delta").template get<double>();
+        // h_min
+        double h_min = constraints.at("h_min").template get<double>();
+        // a
+        double a = constraints.at("a").template get<double>();
+        // b
+        double b = constraints.at("b").template get<double>();
+        // position
+        std::string position =
+            constraints.at("position").template get<std::string>();
+        if (position == "corner")
+          position_ = mpm::Position::Corner;
+        else if (position == "edge")
+          position_ = mpm::Position::Edge;
+        else if (position == "face")
+          position_ = mpm::Position::Face;
+        else
+          position_ = mpm::Position::None;
+        // Add absorbing constraint to mesh
+        auto absorbing_constraint = std::make_shared<mpm::AbsorbingConstraint>(
+            nset_id, dir, delta, h_min, a, b, position_);
+        bool absorbing_constraints =
+            constraints_->assign_nodal_absorbing_constraint(
+                nset_id, absorbing_constraint);
+        if (!absorbing_constraints)
+          throw std::runtime_error(
+              "Nodal absorbing constraint is not properly assigned");
+        absorbing_boundary_ = true;
+        absorbing_nset_id_.emplace_back(nset_id);
+        absorbing_constraint_.emplace_back(absorbing_constraint);
+      }
+    } else
+      throw std::runtime_error("Absorbing constraints JSON not found");
+
+  } catch (std::exception& exception) {
+    console_->warn("#{}: Absorbing conditions are undefined {} ", __LINE__,
+                   exception.what());
+  }
+}
+
+// Apply nodal absorbing constraints
+template <unsigned Tdim>
+void mpm::MPMBase<Tdim>::nodal_absorbing_constraints() {
+  for (int i = 0; i < absorbing_nset_id_.size(); i++) {
+    auto nset_id_ = absorbing_nset_id_.at(i);
+    auto a_constraint_ = absorbing_constraint_.at(i);
+    bool absorbing_constraints =
+        constraints_->assign_nodal_absorbing_constraint(nset_id_,
+                                                        a_constraint_);
+    if (!absorbing_constraints)
+      throw std::runtime_error(
+          "Nodal absorbing constraint is not properly assigned");
+  }
+}
+
 //! Cell entity sets
 template <unsigned Tdim>
 void mpm::MPMBase<Tdim>::cell_entity_sets(const Json& mesh_props,
diff --git a/include/solvers/mpm_explicit.h b/include/solvers/mpm_explicit.h
index 2e463ac2d..1d8abda59 100644
--- a/include/solvers/mpm_explicit.h
+++ b/include/solvers/mpm_explicit.h
@@ -75,6 +75,10 @@ class MPMExplicit : public MPMBase<Tdim> {
   using mpm::MPMBase<Tdim>::damping_factor_;
   //! Locate particles
   using mpm::MPMBase<Tdim>::locate_particles_;
+  //! Constraints Pointer
+  using mpm::MPMBase<Tdim>::constraints_;
+  //! Absorbing Boundary
+  using mpm::MPMBase<Tdim>::absorbing_boundary_;
 
  private:
   //! Pressure smoothing
diff --git a/include/solvers/mpm_explicit.tcc b/include/solvers/mpm_explicit.tcc
index 900a7f27a..bde8799ab 100644
--- a/include/solvers/mpm_explicit.tcc
+++ b/include/solvers/mpm_explicit.tcc
@@ -121,6 +121,7 @@ bool mpm::MPMExplicit<Tdim>::solve() {
   this->initialise_loads();
 
   auto solver_begin = std::chrono::steady_clock::now();
+
   // Main loop
   for (; step_ < nsteps_; ++step_) {
 
@@ -156,6 +157,12 @@ bool mpm::MPMExplicit<Tdim>::solve() {
     mpm_scheme_->compute_forces(gravity_, phase, step_,
                                 set_node_concentrated_force_);
 
+    // Apply Absorbing Constraint
+    if (absorbing_boundary_) {
+      mpm_scheme_->absorbing_boundary_properties();
+      this->nodal_absorbing_constraints();
+    }
+
     // Particle kinematics
     mpm_scheme_->compute_particle_kinematics(velocity_update_, phase, "Cundall",
                                              damping_factor_);
diff --git a/include/solvers/mpm_scheme/mpm_scheme.h b/include/solvers/mpm_scheme/mpm_scheme.h
index 683d149ad..ab08da28d 100644
--- a/include/solvers/mpm_scheme/mpm_scheme.h
+++ b/include/solvers/mpm_scheme/mpm_scheme.h
@@ -57,6 +57,9 @@ class MPMScheme {
       const Eigen::Matrix<double, Tdim, 1>& gravity, unsigned phase,
       unsigned step, bool concentrated_nodal_forces);
 
+  //! Assign relevant properties for absorbing boundary
+  virtual inline void absorbing_boundary_properties();
+
   //! Compute acceleration velocity position
   //! \param[in] velocity_update Velocity or acceleration update flag
   //! \param[in] phase Phase of particle
diff --git a/include/solvers/mpm_scheme/mpm_scheme.tcc b/include/solvers/mpm_scheme/mpm_scheme.tcc
index e95c61795..d558eb4e0 100644
--- a/include/solvers/mpm_scheme/mpm_scheme.tcc
+++ b/include/solvers/mpm_scheme/mpm_scheme.tcc
@@ -168,6 +168,28 @@ inline void mpm::MPMScheme<Tdim>::compute_forces(
 #endif
 }
 
+// Assign Absorbing Boundary Properties
+template <unsigned Tdim>
+inline void mpm::MPMScheme<Tdim>::absorbing_boundary_properties() {
+  mesh_->create_nodal_properties();
+
+  // Initialise nodal properties
+  mesh_->initialise_nodal_properties();
+
+  // Append material ids to nodes
+  mesh_->iterate_over_particles(
+      std::bind(&mpm::ParticleBase<Tdim>::append_material_id_to_nodes,
+                std::placeholders::_1));
+
+  mesh_->iterate_over_particles(
+      std::bind(&mpm::ParticleBase<Tdim>::map_wave_velocities_to_nodes,
+                std::placeholders::_1));
+  // Map multimaterial displacements from particles to nodes
+  mesh_->iterate_over_particles(std::bind(
+      &mpm::ParticleBase<Tdim>::map_multimaterial_displacements_to_nodes,
+      std::placeholders::_1));
+}
+
 // Compute particle kinematics
 template <unsigned Tdim>
 inline void mpm::MPMScheme<Tdim>::compute_particle_kinematics(
diff --git a/tests/interface_test.cc b/tests/interface_test.cc
index cb3b151d7..04ff8b024 100644
--- a/tests/interface_test.cc
+++ b/tests/interface_test.cc
@@ -44,6 +44,9 @@ TEST_CASE("Interface functions are checked", "[interface]") {
                                     Nmaterials);
   nodal_properties->create_property("normal_unit_vectors", Nnodes * Dim,
                                     Nmaterials);
+  nodal_properties->create_property("wave_velocities", Nnodes * Dim,
+                                    Nmaterials);
+  nodal_properties->create_property("density", Nnodes, Nmaterials);
 
   // Element
   std::shared_ptr<mpm::Element<Dim>> element =
@@ -204,6 +207,11 @@ TEST_CASE("Interface functions are checked", "[interface]") {
   node2->compute_multimaterial_normal_unit_vector();
   node3->compute_multimaterial_normal_unit_vector();
 
+  // Map wave velocities to nodes
+  particle1->map_wave_velocities_to_nodes();
+  particle2->map_wave_velocities_to_nodes();
+  particle3->map_wave_velocities_to_nodes();
+
   Eigen::Matrix<double, 4, 2> masses;
   // clang-format off
   masses << 0.96, 1.46,
@@ -283,6 +291,18 @@ TEST_CASE("Interface functions are checked", "[interface]") {
              0.89442719099992,  0.88491822238198;
   // clang-format on
 
+  Eigen::Matrix<double, 8, 2> wave_velocities;
+  // clang-format off
+  wave_velocities << 116.023870223, 116.023870223,
+                     62.0173672946, 62.0173672946,
+                     116.023870223, 116.023870223,
+                     62.0173672946, 62.0173672946,
+                     116.023870223, 116.023870223,
+                     62.0173672946, 62.0173672946,
+                     116.023870223, 116.023870223,
+                     62.0173672946, 62.0173672946;
+  // clang-format on
+
   // Check if nodal properties were properly mapped and computed
   for (int i = 0; i < Nnodes; ++i) {
     for (int j = 0; j < Nmaterials; ++j) {
@@ -304,6 +324,9 @@ TEST_CASE("Interface functions are checked", "[interface]") {
             Approx(gradients(i * Dim + k, j)).epsilon(tolerance));
         REQUIRE(nodal_properties->property("normal_unit_vectors", i, j, Dim)(
                     k, 0) == Approx(normal(i * Dim + k, j)).epsilon(tolerance));
+        REQUIRE(
+            nodal_properties->property("wave_velocities", i, j, Dim)(k, 0) ==
+            Approx(wave_velocities(i * Dim + k, j)).epsilon(tolerance));
       }
       // Check if normal vector are also unit vectors
       REQUIRE(
diff --git a/tests/mesh_test_2d.cc b/tests/mesh_test_2d.cc
index eedc61044..3f93ac785 100644
--- a/tests/mesh_test_2d.cc
+++ b/tests/mesh_test_2d.cc
@@ -1150,6 +1150,93 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") {
                       set_id, friction_constraint) == false);
         }
 
+        // Test assign absorbing constraint
+        SECTION("Check assign absorbing constraints") {
+          tsl::robin_map<mpm::Index, std::vector<mpm::Index>> node_sets;
+          node_sets[0] = std::vector<mpm::Index>{0, 2};
+          node_sets[1] = std::vector<mpm::Index>{1, 3};
+          node_sets[2] = std::vector<mpm::Index>{};
+
+          REQUIRE(mesh->create_node_sets(node_sets, true) == true);
+
+          //! Constraints object
+          auto constraints = std::make_shared<mpm::Constraints<Dim>>(mesh);
+
+          // Assign nodal properties
+          auto nodal_properties = std::make_shared<mpm::NodalProperties>();
+          nodal_properties->create_property("density", 1, 1);
+          nodal_properties->create_property("wave_velocities", Dim, 1);
+          nodal_properties->create_property("displacements", Dim, 1);
+          for (double i = 0; i < 4; ++i) {
+            REQUIRE_NOTHROW(
+                mesh->node(i)->initialise_property_handle(0, nodal_properties));
+            mesh->node(i)->append_material_id(0);
+          }
+
+          int set_id = 0;
+          unsigned dir = 0;
+          double delta = 1;
+          double h_min = 1;
+          double a = 1;
+          double b = 1;
+          mpm::Position pos = mpm::Position::Edge;
+
+          // Add absorbing constraint to mesh
+          auto absorbing_constraint =
+              std::make_shared<mpm::AbsorbingConstraint>(set_id, dir, delta,
+                                                         h_min, a, b, pos);
+          REQUIRE(constraints->assign_nodal_absorbing_constraint(
+                      set_id, absorbing_constraint) == true);
+
+          set_id = 1;
+          dir = 1;
+          delta = 3;
+          h_min = 12;
+          a = 2;
+          b = 2;
+          pos = mpm::Position::Corner;
+          // Add absorbing constraint to mesh
+          absorbing_constraint = std::make_shared<mpm::AbsorbingConstraint>(
+              set_id, dir, delta, h_min, a, b, pos);
+          REQUIRE(constraints->assign_nodal_absorbing_constraint(
+                      set_id, absorbing_constraint) == true);
+
+          // When constraints fail: invalid direction
+          dir = 3;
+          // Add absorbing constraint to mesh
+          absorbing_constraint = std::make_shared<mpm::AbsorbingConstraint>(
+              set_id, dir, delta, h_min, a, b, pos);
+          REQUIRE(constraints->assign_nodal_absorbing_constraint(
+                      set_id, absorbing_constraint) == false);
+
+          // When constraints fail: invalid delta
+          dir = 1;
+          delta = 0;
+          // Add absorbing constraint to mesh
+          absorbing_constraint = std::make_shared<mpm::AbsorbingConstraint>(
+              set_id, dir, delta, h_min, a, b, pos);
+          REQUIRE(constraints->assign_nodal_absorbing_constraint(
+                      set_id, absorbing_constraint) == false);
+
+          // When constraints fail: empty node set
+          delta = 3;
+          set_id = 2;
+          // Add absorbing constraint to mesh
+          absorbing_constraint = std::make_shared<mpm::AbsorbingConstraint>(
+              set_id, dir, delta, h_min, a, b, pos);
+          REQUIRE(constraints->assign_nodal_absorbing_constraint(
+                      set_id, absorbing_constraint) == false);
+
+          // When constraints fail: invalid absorbing boundary position
+          set_id = 1;
+          pos = mpm::Position::None;
+          // Add absorbing constraint to mesh
+          absorbing_constraint = std::make_shared<mpm::AbsorbingConstraint>(
+              set_id, dir, delta, h_min, a, b, pos);
+          REQUIRE(constraints->assign_nodal_absorbing_constraint(
+                      set_id, absorbing_constraint) == false);
+        }
+
         // Test assign velocity constraints to nodes
         SECTION("Check assign velocity constraints to nodes") {
           // Vector of particle coordinates
@@ -1192,6 +1279,60 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") {
                       friction_constraints) == false);
         }
 
+        // Test assign absorbing constraints to nodes
+        SECTION("Check assign absorbing constraints to nodes") {
+          // Vector of particle coordinates
+          std::vector<std::tuple<mpm::Index, unsigned, double, double, double,
+                                 double, mpm::Position>>
+              absorbing_constraints;
+          //! Constraints object
+          auto constraints = std::make_shared<mpm::Constraints<Dim>>(mesh);
+
+          // Assign nodal properties
+          auto nodal_properties = std::make_shared<mpm::NodalProperties>();
+          nodal_properties->create_property("density", 1, 1);
+          nodal_properties->create_property("wave_velocities", Dim, 1);
+          nodal_properties->create_property("displacements", Dim, 1);
+          for (double i = 0; i < 4; ++i) {
+            REQUIRE_NOTHROW(
+                mesh->node(i)->initialise_property_handle(0, nodal_properties));
+          }
+
+          // Constraint
+          absorbing_constraints.emplace_back(
+              std::make_tuple(0, 0, 1, 3, 2, 2, mpm::Position::Edge));
+          absorbing_constraints.emplace_back(
+              std::make_tuple(1, 1, 2, 4, 1, 1, mpm::Position::Edge));
+          absorbing_constraints.emplace_back(
+              std::make_tuple(2, 0, 1, 1, 1, 1, mpm::Position::Edge));
+          absorbing_constraints.emplace_back(
+              std::make_tuple(3, 1, 3, 2, 3, 3, mpm::Position::Edge));
+
+          REQUIRE(constraints->assign_nodal_absorbing_constraints(
+                      absorbing_constraints) == true);
+
+          // When constraints fail: invalid direction
+          absorbing_constraints.clear();
+          absorbing_constraints.emplace_back(
+              std::make_tuple(3, 2, 3, 2, 3, 3, mpm::Position::Edge));
+          REQUIRE(constraints->assign_nodal_absorbing_constraints(
+                      absorbing_constraints) == false);
+
+          // When constraints fail: invalid delta
+          absorbing_constraints.clear();
+          absorbing_constraints.emplace_back(
+              std::make_tuple(3, 1, 1, 3, 1, 1, mpm::Position::Edge));
+          REQUIRE(constraints->assign_nodal_absorbing_constraints(
+                      absorbing_constraints) == false);
+
+          // When constraints fail: invalid position
+          absorbing_constraints.clear();
+          absorbing_constraints.emplace_back(
+              std::make_tuple(0, 0, 1, 3, 2, 2, mpm::Position::None));
+          REQUIRE(constraints->assign_nodal_absorbing_constraints(
+                      absorbing_constraints) == false);
+        }
+
         // Test assign nodes concentrated_forces
         SECTION("Check assign nodes concentrated_forces") {
           // Vector of node coordinates
diff --git a/tests/mesh_test_3d.cc b/tests/mesh_test_3d.cc
index 06f4e0cf2..e6698f764 100644
--- a/tests/mesh_test_3d.cc
+++ b/tests/mesh_test_3d.cc
@@ -1261,6 +1261,67 @@ TEST_CASE("Mesh is checked for 3D case", "[mesh][3D]") {
                       set_id, friction_constraint) == false);
         }
 
+        // Test assign absorbing constraint
+        SECTION("Check assign absorbing constraints") {
+          tsl::robin_map<mpm::Index, std::vector<mpm::Index>> node_sets;
+          node_sets[0] = std::vector<mpm::Index>{0, 2};
+          node_sets[1] = std::vector<mpm::Index>{1, 3};
+          node_sets[2] = std::vector<mpm::Index>{};
+
+          REQUIRE(mesh->create_node_sets(node_sets, true) == true);
+
+          //! Constraints object
+          auto constraints = std::make_shared<mpm::Constraints<Dim>>(mesh);
+
+          // Assign nodal properties
+          auto nodal_properties = std::make_shared<mpm::NodalProperties>();
+          nodal_properties->create_property("density", 1, 1);
+          nodal_properties->create_property("wave_velocities", Dim, 1);
+          nodal_properties->create_property("displacements", Dim, 1);
+          for (double i = 0; i < 4; ++i) {
+            REQUIRE_NOTHROW(
+                mesh->node(i)->initialise_property_handle(0, nodal_properties));
+            mesh->node(i)->append_material_id(0);
+          }
+
+          int set_id = 0;
+          unsigned dir = 0;
+          double delta = 1;
+          double h_min = 1;
+          double a = 1;
+          double b = 1;
+          mpm::Position pos = mpm::Position::Edge;
+
+          // Add absorbing constraint to mesh
+          auto absorbing_constraint =
+              std::make_shared<mpm::AbsorbingConstraint>(set_id, dir, delta,
+                                                         h_min, a, b, pos);
+          REQUIRE(constraints->assign_nodal_absorbing_constraint(
+                      set_id, absorbing_constraint) == true);
+
+          pos = mpm::Position::Corner;
+          // Add absorbing constraint to mesh
+          absorbing_constraint = std::make_shared<mpm::AbsorbingConstraint>(
+              set_id, dir, delta, h_min, a, b, pos);
+          REQUIRE(constraints->assign_nodal_absorbing_constraint(
+                      set_id, absorbing_constraint) == true);
+
+          pos = mpm::Position::Face;
+          // Add absorbing constraint to mesh
+          absorbing_constraint = std::make_shared<mpm::AbsorbingConstraint>(
+              set_id, dir, delta, h_min, a, b, pos);
+          REQUIRE(constraints->assign_nodal_absorbing_constraint(
+                      set_id, absorbing_constraint) == true);
+
+          // When constraints fail: invalid absorbing boundary position
+          pos = mpm::Position::None;
+          // Add absorbing constraint to mesh
+          absorbing_constraint = std::make_shared<mpm::AbsorbingConstraint>(
+              set_id, dir, delta, h_min, a, b, pos);
+          REQUIRE(constraints->assign_nodal_absorbing_constraint(
+                      set_id, absorbing_constraint) == false);
+        }
+
         // Test assign velocity constraints to nodes
         SECTION("Check assign velocity constraints to nodes") {
           // Vector of particle coordinates