diff --git a/examples/acoustic_2D.py b/examples/acoustic_2D.py
index 9641f47..5a9f67d 100644
--- a/examples/acoustic_2D.py
+++ b/examples/acoustic_2D.py
@@ -61,19 +61,19 @@
 # config boundary conditions
 # (none,  null_dirichlet or null_neumann)
 space_model.config_boundary(
-    damping_length=0,
+    damping_length=(0, 510, 510, 510),
     boundary_condition=(
         "null_neumann", "null_dirichlet",
-        "none", "null_dirichlet"
-    ),
-    damping_polynomial_degree=3,
-    damping_alpha=0.001
+        "null_dirichlet", "null_dirichlet"
+    )
 )
 
+print(' damping_alpha=',space_model.damping_alpha)
+
 # create the time model
 time_model = TimeModel(
     space_model=space_model,
-    tf=1.0,
+    tf=2.0,
     saving_stride=0
 )
 
diff --git a/simwave/kernel/frontend/model.py b/simwave/kernel/frontend/model.py
index 91fdb65..7075472 100644
--- a/simwave/kernel/frontend/model.py
+++ b/simwave/kernel/frontend/model.py
@@ -187,7 +187,24 @@ def damping_alpha(self):
         try:
             return self._damping_alpha
         except AttributeError:
-            return 0.001
+            if self.dimension == 2:
+                # 2 dimension
+                z_min, z_max, x_min, x_max = self.bounding_box
+                z_top_damp, z_bottom_damp, x_left_damp, x_right_damp = self.damping_length
+                nz = z_max - z_min + z_top_damp + z_bottom_damp
+                nx = x_max - x_min + x_left_damp + x_right_damp
+                # get the maximum domain length in one of 2 dimensions (meters)
+                max_length = max(nz, nx)
+            else:
+                # 3 dimension
+                z_min, z_max, x_min, x_max, y_min, y_max = self.bounding_box
+                z_top_damp, z_bottom_damp, x_left_damp, x_right_damp, y_front_damp, y_back_damp = self.damping_length
+                nz = z_max - z_min + z_top_damp + z_bottom_damp
+                nx = x_max - x_min + x_left_damp + x_right_damp
+                ny = y_max - y_min + y_front_damp + y_back_damp
+                # get the maximum domain length in one of 3 dimensions (meters)
+                max_length = max(nz, nx, ny)
+            return np.log(1/10 ** -3) * (self.dimension * np.max(self.velocity_model) / (2 * max_length) ) / max_length ** 2
 
     def fd_coefficients(self, derivative_order):
         """
@@ -269,7 +286,7 @@ def halo_size(self):
         return (space_order_radius, ) * self.dimension * 2
 
     def config_boundary(self, damping_length=0.0, boundary_condition="none",
-                        damping_polynomial_degree=3, damping_alpha=0.001):
+                        damping_polynomial_degree=2, damping_alpha=None):
         """
         Applies the domain extension (for absorbing layers with damping)
         and boundary conditions.
@@ -296,7 +313,8 @@ def config_boundary(self, damping_length=0.0, boundary_condition="none",
         """
 
         self._damping_polynomial_degree = damping_polynomial_degree
-        self._damping_alpha = damping_alpha
+        if damping_alpha is not None: 
+            self._damping_alpha = damping_alpha
 
         # if it is not, convert damping_length to tuple
         if isinstance(damping_length, (float, int)):
@@ -394,11 +412,9 @@ def damping_mask(self):
             mode="linear_ramp",
             end_values=self.nbl_pad_width
         )
-
         # change the damping values (coefficients) according to a function
         degree = self.damping_polynomial_degree
         damp_mask = (damp_mask ** degree) * self.damping_alpha
-
         # damp mask in the halo region
         # The values in this extended region is zero
         damp_mask = np.pad(array=damp_mask, pad_width=self.halo_pad_width)