diff --git a/CHANGELOG.md b/CHANGELOG.md
index b3a5bf8ff..9b59c83d6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -18,6 +18,7 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm
 - The World Builder Visualizer will now use zlib compression for vtu files by default. If zlib is not available binary output will be used. \[Menno Fraters; 2021-06-26; [#282](github.com/GeodynamicWorldBuilder/WorldBuilder/pull/282)\]
 - The return argument type of the distance_point_from_curved_planes function has been converted from a map to a struct, requiring a change in the plugin interfaces for 'fault_models' and 'subducting_plate_models', but significantly increasing the speed of the function and all functions that access its returned values. \[Rene Gassmoeller; 2021-06-27; [#289](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/289)\]
 - The plugin systems (temperature, composition and grains) and the distance_point_from_curved_planes function now all pass a precomputed NaturalCoordinate, besides just the cartesian position. It turns out that this can make a significant performance differce. \[Issue found and solution suggested by Wolfgang Bangerth, implemented and tested by Menno Fraters; 2021-07-03; [#300](github.com/GeodynamicWorldBuilder/WorldBuilder/pull/300) and [#219](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/219)\]
+- Introduces a bounding box for searching the fault and the subducting plate model properties (temperature, composition, grains) in the lateral direction. This reduces the computation time by constraining the number of points for which the expensive distance_point_from_curved_planes function is called. \[Arushi Saxena; 2021-07-15; [#219](https://github.com/GeodynamicWorldBuilder/WorldBuilder/issues/219)\]
 
 ## [0.4.0]
 ### Added
diff --git a/include/world_builder/features/fault.h b/include/world_builder/features/fault.h
index 9d0c28641..8e86633f2 100644
--- a/include/world_builder/features/fault.h
+++ b/include/world_builder/features/fault.h
@@ -25,6 +25,7 @@
 #include "world_builder/features/fault_models/grains/interface.h"
 #include "world_builder/features/fault_models/temperature/interface.h"
 #include "world_builder/types/segment.h"
+#include "world_builder/bounding_box.h"
 
 
 namespace WorldBuilder
@@ -104,6 +105,18 @@ namespace WorldBuilder
                            const unsigned int composition_number,
                            double composition_value) const override final;
 
+        /**
+         * Computes the bounding points for a BoundingBox object using two extreme points in all the surface
+         * coordinates and an additional buffer zone that accounts for the fault thickness and length. The first and second
+         * points correspond to the lower left and the upper right corners of the bounding box, respectively (see the
+         * documentation in include/bounding_box.h).
+         * For the spherical system, the buffer zone along the longitudal direction is calculated using the
+         * correponding latitude points.
+         */
+
+        BoundingBox<2>  get_bounding_box (const WorldBuilder::Utilities::NaturalCoordinate &position_in_natural_coordinates,
+                                          const double depth) const;
+
         /**
          * Returns a grains (rotation matrix and grain size) based on the
          * given position, depth in the model, the composition (e.g. representing
@@ -179,6 +192,13 @@ namespace WorldBuilder
         double maximum_total_fault_length;
         double maximum_fault_thickness;
 
+        double min_along_x;
+        double max_along_x;
+        double min_along_y;
+        double max_along_y;
+        double min_lat_cos_inv;
+        double max_lat_cos_inv;
+        double buffer_around_fault_cartesian;
 
     };
   } // namespace Features
diff --git a/include/world_builder/features/subducting_plate.h b/include/world_builder/features/subducting_plate.h
index bf3e9d74d..8f983ce69 100644
--- a/include/world_builder/features/subducting_plate.h
+++ b/include/world_builder/features/subducting_plate.h
@@ -25,6 +25,7 @@
 #include "world_builder/features/subducting_plate_models/grains/interface.h"
 #include "world_builder/features/subducting_plate_models/temperature/interface.h"
 #include "world_builder/types/segment.h"
+#include "world_builder/bounding_box.h"
 
 
 namespace WorldBuilder
@@ -83,6 +84,17 @@ namespace WorldBuilder
         void parse_entries(Parameters &prm) override final;
 
 
+        /**
+         * Computes the bounding points for a BoundingBox object using two extreme points in all the surface
+         * coordinates and an additional buffer zone that accounts for the fault thickness and length. The first and second
+         * points correspond to the lower left and the upper right corners of the bounding box, respectively (see the
+         * documentation in include/bounding_box.h).
+         * For the spherical system, the buffer zone along the longitudal direction is calculated using the
+         * correponding latitude points.
+         */
+        BoundingBox<2>  get_bounding_box (const WorldBuilder::Utilities::NaturalCoordinate &position_in_natural_coordinates,
+                                          const double depth) const;
+
 
         /**
          * Returns a temperature based on the given position, depth in the model,
@@ -182,6 +194,13 @@ namespace WorldBuilder
         double maximum_total_slab_length;
         double maximum_slab_thickness;
 
+        double min_along_x;
+        double max_along_x;
+        double min_along_y;
+        double max_along_y;
+        double min_lat_cos_inv;
+        double max_lat_cos_inv;
+        double buffer_around_slab_cartesian;
     };
   } // namespace Features
 } // namespace WorldBuilder
diff --git a/source/features/fault.cc b/source/features/fault.cc
index 7937d9659..02aaadd91 100644
--- a/source/features/fault.cc
+++ b/source/features/fault.cc
@@ -26,6 +26,7 @@
 #include "world_builder/types/point.h"
 #include "world_builder/types/unsigned_int.h"
 #include "world_builder/world.h"
+#include <algorithm>
 
 using namespace std;
 
@@ -325,6 +326,59 @@ namespace WorldBuilder
           total_fault_length[i] = local_total_fault_length;
           maximum_total_fault_length = std::max(maximum_total_fault_length, local_total_fault_length);
         }
+
+
+      std::vector<double> x_list(original_number_of_coordinates,0.0);
+      std::vector<double> y_list(original_number_of_coordinates,0.0);
+
+      for (size_t j=0; j<original_number_of_coordinates; ++j)
+        {
+          x_list[j] = coordinates[j][0];
+          y_list[j] = coordinates[j][1];
+        }
+
+
+      min_along_x = *std::min_element(x_list.begin(), x_list.end());
+      max_along_x = *std::max_element(x_list.begin(), x_list.end());
+      min_along_y = *std::min_element(y_list.begin(), y_list.end());
+      max_along_y = *std::max_element(y_list.begin(), y_list.end());
+
+      min_lat_cos_inv = 1 / std::cos(min_along_y);
+      max_lat_cos_inv = 1 / std::cos(max_along_y);
+
+      buffer_around_fault_cartesian =  (maximum_fault_thickness + maximum_total_fault_length);
+
+    }
+
+
+    BoundingBox<2>
+    Fault::get_bounding_box (const NaturalCoordinate &position_in_natural_coordinates,
+                             const double depth) const
+    {
+      BoundingBox<2> surface_bounding_box;
+      const double starting_radius_inv = 1 / (position_in_natural_coordinates.get_depth_coordinate() + depth - starting_depth);
+      if (world->parameters.coordinate_system->natural_coordinate_system() == CoordinateSystem::spherical)
+        {
+          std::pair<Point<2>, Point<2> > &spherical_bounding_box = surface_bounding_box.get_boundary_points();
+
+          const double buffer_around_fault_spherical = 2 * const_pi * buffer_around_fault_cartesian * starting_radius_inv;
+
+          spherical_bounding_box.first = {(min_along_x - buffer_around_fault_spherical * min_lat_cos_inv) ,
+                                          (min_along_y - buffer_around_fault_spherical), spherical
+                                         } ;
+
+          spherical_bounding_box.second = {(max_along_x + buffer_around_fault_spherical * max_lat_cos_inv) ,
+                                           (max_along_y + buffer_around_fault_spherical), spherical
+                                          };
+        }
+      else if (world->parameters.coordinate_system->natural_coordinate_system() == CoordinateSystem::cartesian)
+        {
+          std::pair<Point<2>, Point<2> > &bounding_box = surface_bounding_box.get_boundary_points();
+          bounding_box.first = {min_along_x, min_along_y, cartesian};
+          bounding_box.second = {max_along_x, max_along_x, cartesian};
+          surface_bounding_box.extend(buffer_around_fault_cartesian);
+        }
+      return surface_bounding_box;
     }
 
 
@@ -348,7 +402,9 @@ namespace WorldBuilder
               );
 
       // todo: explain and check -starting_depth
-      if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_fault_length + maximum_fault_thickness)
+      if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_fault_length + maximum_fault_thickness &&
+          get_bounding_box(position_in_natural_coordinates, depth).point_inside(Point<2>(position_in_natural_coordinates.get_surface_coordinates(),
+                                                                                world->parameters.coordinate_system->natural_coordinate_system())))
         {
           // todo: explain
           // This function only returns positive values, because we want
@@ -468,6 +524,7 @@ namespace WorldBuilder
       return temperature;
     }
 
+
     double
     Fault::composition(const Point<3> &position_in_cartesian_coordinates,
                        const NaturalCoordinate &position_in_natural_coordinates,
@@ -479,7 +536,9 @@ namespace WorldBuilder
       const double starting_radius = position_in_natural_coordinates.get_depth_coordinate() + depth - starting_depth;
 
       // todo: explain and check -starting_depth
-      if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_fault_length + maximum_fault_thickness)
+      if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_fault_length + maximum_fault_thickness &&
+          get_bounding_box(position_in_natural_coordinates, depth).point_inside(Point<2>(position_in_natural_coordinates.get_surface_coordinates(),
+                                                                                world->parameters.coordinate_system->natural_coordinate_system())))
         {
           // todo: explain
           // This function only returns positive values, because we want
@@ -613,7 +672,9 @@ namespace WorldBuilder
       const double starting_radius = position_in_natural_coordinates.get_depth_coordinate() + depth - starting_depth;
 
       // todo: explain and check -starting_depth
-      if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_fault_length + maximum_fault_thickness)
+      if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_fault_length + maximum_fault_thickness &&
+          get_bounding_box(position_in_natural_coordinates, depth).point_inside(Point<2>(position_in_natural_coordinates.get_surface_coordinates(),
+                                                                                world->parameters.coordinate_system->natural_coordinate_system())))
         {
           // todo: explain
           // This function only returns positive values, because we want
diff --git a/source/features/subducting_plate.cc b/source/features/subducting_plate.cc
index 67be1c3a7..971e40eaa 100644
--- a/source/features/subducting_plate.cc
+++ b/source/features/subducting_plate.cc
@@ -26,6 +26,7 @@
 #include "world_builder/types/point.h"
 #include "world_builder/types/unsigned_int.h"
 #include "world_builder/world.h"
+#include <algorithm>
 
 using namespace std;
 
@@ -329,6 +330,62 @@ namespace WorldBuilder
           total_slab_length[i] = local_total_slab_length;
           maximum_total_slab_length = std::max(maximum_total_slab_length, local_total_slab_length);
         }
+
+      std::vector<double> x_list(original_number_of_coordinates,0.0);
+      std::vector<double> y_list(original_number_of_coordinates,0.0);
+
+      for (size_t j=0; j<original_number_of_coordinates; ++j)
+        {
+          x_list[j] = coordinates[j][0];
+          y_list[j] = coordinates[j][1];
+        }
+
+      // Here, we compute the spherical bounding box using the two extreme points of the box containing all the surface
+      // coordinates and an additional buffer zone that accounts for the fault thickness and length. The first and second
+      // points correspond to the lower left and the upper right corners of the bounding box, respectively (see the
+      // documentation in include/bounding_box.h).
+      // For the spherical system, the buffer zone along the longitudal direction is calculated using the
+      // correponding latitude points.
+
+      min_along_x = *std::min_element(x_list.begin(), x_list.end());
+      max_along_x = *std::max_element(x_list.begin(), x_list.end());
+      min_along_y = *std::min_element(y_list.begin(), y_list.end());
+      max_along_y = *std::max_element(y_list.begin(), y_list.end());
+
+      min_lat_cos_inv = 1 / std::cos(min_along_y);
+      max_lat_cos_inv = 1 / std::cos(max_along_y);
+
+      buffer_around_slab_cartesian =  (maximum_slab_thickness + maximum_total_slab_length);
+    }
+
+    BoundingBox<2>
+    SubductingPlate::get_bounding_box (const NaturalCoordinate &position_in_natural_coordinates,
+                                       const double depth) const
+    {
+      BoundingBox<2> surface_bounding_box;
+      const double starting_radius_inv = 1 / (position_in_natural_coordinates.get_depth_coordinate() + depth - starting_depth);
+      if (world->parameters.coordinate_system->natural_coordinate_system() == CoordinateSystem::spherical)
+        {
+          std::pair<Point<2>, Point<2> > &spherical_bounding_box = surface_bounding_box.get_boundary_points();
+
+          const double buffer_around_fault_spherical = 2 * const_pi * buffer_around_slab_cartesian * starting_radius_inv;
+
+          spherical_bounding_box.first = {(min_along_x - buffer_around_fault_spherical * min_lat_cos_inv) ,
+                                          (min_along_y - buffer_around_fault_spherical), spherical
+                                         } ;
+
+          spherical_bounding_box.second = {(max_along_x + buffer_around_fault_spherical * max_lat_cos_inv) ,
+                                           (max_along_y + buffer_around_fault_spherical), spherical
+                                          };
+        }
+      else if (world->parameters.coordinate_system->natural_coordinate_system() == CoordinateSystem::cartesian)
+        {
+          std::pair<Point<2>, Point<2> > &bounding_box = surface_bounding_box.get_boundary_points();
+          bounding_box.first = {min_along_x, min_along_y, cartesian};
+          bounding_box.second = {max_along_x, max_along_x, cartesian};
+          surface_bounding_box.extend(buffer_around_slab_cartesian);
+        }
+      return surface_bounding_box;
     }
 
 
@@ -352,7 +409,9 @@ namespace WorldBuilder
               );
 
       // todo: explain and check -starting_depth
-      if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_slab_length + maximum_slab_thickness)
+      if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_slab_length + maximum_slab_thickness &&
+          get_bounding_box(position_in_natural_coordinates, depth).point_inside(Point<2>(position_in_natural_coordinates.get_surface_coordinates(),
+                                                                                world->parameters.coordinate_system->natural_coordinate_system())))
         {
           /*WBAssert(coordinates.size() == slab_segment_lengths.size(),
                    "Internal error: The size of coordinates (" << coordinates.size()
@@ -489,7 +548,9 @@ namespace WorldBuilder
       const double starting_radius = position_in_natural_coordinates.get_depth_coordinate() + depth - starting_depth;
 
       // todo: explain and check -starting_depth
-      if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_slab_length + maximum_slab_thickness)
+      if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_slab_length + maximum_slab_thickness &&
+          get_bounding_box(position_in_natural_coordinates, depth).point_inside(Point<2>(position_in_natural_coordinates.get_surface_coordinates(),
+                                                                                world->parameters.coordinate_system->natural_coordinate_system())))
         {
           // todo: explain
           WorldBuilder::Utilities::PointDistanceFromCurvedPlanes distance_from_planes =
@@ -620,7 +681,9 @@ namespace WorldBuilder
       const double starting_radius = position_in_natural_coordinates.get_depth_coordinate() + depth - starting_depth;
 
       // todo: explain and check -starting_depth
-      if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_slab_length + maximum_slab_thickness)
+      if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_slab_length + maximum_slab_thickness &&
+          get_bounding_box(position_in_natural_coordinates, depth).point_inside(Point<2>(position_in_natural_coordinates.get_surface_coordinates(),
+                                                                                world->parameters.coordinate_system->natural_coordinate_system())))
         {
           // todo: explain
           WorldBuilder::Utilities::PointDistanceFromCurvedPlanes distance_from_planes =