diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md
index 2729772b8..d8c279328 100755
--- a/docs/CHANGELOG.md
+++ b/docs/CHANGELOG.md
@@ -1,6 +1,16 @@
 All notable changes to this project will be documented in this file.
 We follow the [Semantic Versioning 2.0.0](http://semver.org/) format.
 
+## v4.8.10.7 - 2025-09-19 - [PR#1608](https://github.com/NOAA-OWP/inundation-mapping/pull/1608)
+
+Changes files to use `fiona` as the I/O engine when reading files with `gpd.read_file()` to avoid segmentation faults. The issue was originally documented in #1376 and fixed in #1490, but two more files are fixed here.
+
+### Changes
+
+- `src/`
+    - `mask_dem.py` and `mitigate_branch_outlet_backpool.py`: Added `engine='fiona'` to `gpd.read_file()`. Also removed an unnecessary file read in `mask_dem.py`.
+<br/>
+
 ## v4.8.10.6 - 2025-09-19 - [PR#1642](https://github.com/NOAA-OWP/inundation-mapping/pull/1642)
 
 Adjust eval_plot.py so it only looks for the config file when in -sp (spatial) mode. Spatial mode is used when running the tool to create FIM_performance files. In that mode, it can only run in OWP servers as it needs to talk to internal servers.
@@ -279,10 +289,9 @@ Addresses bug related to the `location_id` data type that is read in from the `a
 
 ### Changes
 `src/src_adjust_usgs_rating_trace.py`: Added `dtype={'location_id': object}` to the acceptable_sites csv file read
-  
+ 
 <br/><br/>
 
-
 ## v4.8.7.1 - 2025-07-18 - [PR#1539]([https://github.com/NOAA-OWP/inundation-mapping/pull/1539])
 
 Updates the CatFIM site comparison tool to make the outputs better suited to be loaded into HydroVis. Add % change calculations to site change outputs. 
diff --git a/src/mask_dem.py b/src/mask_dem.py
index 7f92ffc7a..c1ea77ece 100755
--- a/src/mask_dem.py
+++ b/src/mask_dem.py
@@ -11,7 +11,7 @@
 from rasterio.mask import mask
 
 
-gpd.options.io_engine = "pyogrio"
+# gpd.options.io_engine = "pyogrio"
 
 
 def mask_dem(
@@ -71,7 +71,9 @@ def mask_dem(
 
         elif os.path.exists(levee_levelpaths):
             # Mask levee-protected areas protected against level path
+            catchments = gpd.read_file(catchments_filename, engine='fiona')
             levee_levelpaths = pd.read_csv(levee_levelpaths)
+            leveed = gpd.read_file(nld_filename, engine='fiona')
 
             # Select levees associated with branch
             levee_levelpaths = levee_levelpaths[levee_levelpaths[branch_id_attribute] == branch_id]
@@ -80,8 +82,6 @@ def mask_dem(
             levelpath_levees = list(levee_levelpaths[levee_id_attribute])
 
             if len(levelpath_levees) > 0:
-                leveed = gpd.read_file(nld_filename, engine='fiona')
-
                 # Get geometries of levee protected areas associated with levelpath
                 geoms = [
                     feature['geometry']
@@ -93,9 +93,6 @@ def mask_dem(
                     dem_masked, _ = mask(dem, geoms, invert=True)
 
             # Mask levee-protected areas not protected against level path
-            catchments = gpd.read_file(catchments_filename, engine='fiona')
-            leveed = gpd.read_file(nld_filename, engine='fiona')
-
             leveed_area_catchments = gpd.overlay(catchments, leveed, how="union")
 
             # Select levee catchments not associated with level path
diff --git a/src/mitigate_branch_outlet_backpool.py b/src/mitigate_branch_outlet_backpool.py
index 607fc57ab..e99014f10 100755
--- a/src/mitigate_branch_outlet_backpool.py
+++ b/src/mitigate_branch_outlet_backpool.py
@@ -266,7 +266,7 @@ def mask_raster_to_boundary(raster_path, boundary_json, save_path):
 
     # --------------------------------------------------------------
     # Read in nwm lines, explode to ensure linestrings are the only geometry
-    nwm_streams = gpd.read_file(nwm_streams_filename).explode(index_parts=True)
+    nwm_streams = gpd.read_file(nwm_streams_filename, engine='fiona').explode(index_parts=True)
 
     # Check whether it's branch zero
     if 'levpa_id' in nwm_streams.columns:
@@ -294,8 +294,8 @@ def mask_raster_to_boundary(raster_path, boundary_json, save_path):
         #     print(f'No catchment pixels geome at {catchment_reaches_filename}.')
 
         # Read in split_flows_file and split_points_filename
-        split_flows_geom = gpd.read_file(split_flows_filename)
-        split_points_geom = gpd.read_file(split_points_filename)
+        split_flows_geom = gpd.read_file(split_flows_filename, engine='fiona')
+        split_points_geom = gpd.read_file(split_points_filename, engine='fiona')
 
         # Subset the split flows to get the last one
         split_flows_last_geom = split_flows_geom[split_flows_geom['NextDownID'] == '-1']
@@ -423,7 +423,7 @@ def mask_raster_to_boundary(raster_path, boundary_json, save_path):
                     # print("gdal_polygonize executed successfully.")  # verbose
 
                     # Read in the polygonized catchment pixels
-                    cp_poly_geom = gpd.read_file(catchment_pixels_polygonized_filename)
+                    cp_poly_geom = gpd.read_file(catchment_pixels_polygonized_filename, engine='fiona')
 
                     # --------------------------------------------------------------
                     # Mask problematic pixel catchment from the catchments rasters