diff --git a/stardis/radiation_field/opacities/opacities_solvers/util.py b/stardis/radiation_field/opacities/opacities_solvers/util.py
index 7e297e80..0d83b2cb 100644
--- a/stardis/radiation_field/opacities/opacities_solvers/util.py
+++ b/stardis/radiation_field/opacities/opacities_solvers/util.py
@@ -1,5 +1,6 @@
 import numpy as np
 import pandas as pd
+import logging
 
 from astropy import units as u, constants as const
 from tardis.util.base import species_string_to_tuple
@@ -7,6 +8,9 @@
 from scipy.interpolate import LinearNDInterpolator
 
 
+logger = logging.getLogger(__name__)
+
+
 def sigma_file(tracing_lambdas, temperatures, fpath, opacity_source=None):
     """
     Reads and interpolates a cross-section file.
@@ -46,12 +50,16 @@ def sigma_file(tracing_lambdas, temperatures, fpath, opacity_source=None):
         linear_interp_h2plus_bf = LinearNDInterpolator(
             np.vstack([file_waves_mesh.ravel(), file_temps_mesh.ravel()]).T,
             file_cross_sections.flatten(),
+            fill_value=0,
         )
         lambdas, temps = np.meshgrid(tracing_lambdas, temperatures)
         sigmas = (
             linear_interp_h2plus_bf(lambdas, temps) * 1e-18
         )  # Scaling from Stancil 1994 table
-
+        if np.any(sigmas == 0):
+            logger.warning(
+                f"Outside of interpolation range for H2+ BF cross-sections at depth points {np.unique(np.where(sigmas == 0)[0])}. Assuming 0 opacity from H-FF for these depth points."
+            )
     elif (
         opacity_source == "Hminus_ff"
     ):  # This section specifically ingests the Bell and Berrington 1987 h_minus_ff_B1987.dat table found in data.
@@ -68,6 +76,7 @@ def sigma_file(tracing_lambdas, temperatures, fpath, opacity_source=None):
         linear_interp_hminus_ff = LinearNDInterpolator(
             np.vstack([file_waves_mesh.ravel(), file_thetas_mesh.ravel()]).T,
             file_values.flatten(),
+            fill_value=0,
         )
         lambdas, thetas = np.meshgrid(tracing_lambdas, 5040 / temperatures)
         sigmas = (
@@ -76,6 +85,10 @@ def sigma_file(tracing_lambdas, temperatures, fpath, opacity_source=None):
             * const.k_B.cgs.value
             * temperatures[:, np.newaxis]
         )
+        if np.any(sigmas == 0):
+            logger.warning(
+                f"Outside of interpolation range for H- FF cross-sections at depth points {np.unique(np.where(sigmas == 0)[0])}. Assuming 0 opacity from H-FF for these depth points."
+            )
 
     elif (
         opacity_source == "Hminus_bf"