diff --git a/crates/bevy_pbr/src/atmosphere/aerial_view_lut.wgsl b/crates/bevy_pbr/src/atmosphere/aerial_view_lut.wgsl
index 4f9a05e56a360..cb37d774b2b27 100644
--- a/crates/bevy_pbr/src/atmosphere/aerial_view_lut.wgsl
+++ b/crates/bevy_pbr/src/atmosphere/aerial_view_lut.wgsl
@@ -7,7 +7,7 @@
             sample_transmittance_lut, sample_density_lut, rayleigh, henyey_greenstein,
             sample_multiscattering_lut, AtmosphereSample, sample_local_inscattering,
             uv_to_ndc, max_atmosphere_distance, uv_to_ray_direction, 
-            MIDPOINT_RATIO, get_view_position
+            MIDPOINT_RATIO, get_view_position, MIN_EXTINCTION
         },
     }
 }
@@ -52,7 +52,7 @@ fn main(@builtin(global_invocation_id) idx: vec3<u32>) {
             var inscattering = sample_local_inscattering(scattering, ray_dir.xyz, sample_pos);
 
             // Analytical integration of the single scattering term in the radiance transfer equation
-            let s_int = (inscattering - inscattering * sample_transmittance) / extinction;
+            let s_int = (inscattering - inscattering * sample_transmittance) / max(extinction, MIN_EXTINCTION);
             total_inscattering += throughput * s_int;
 
             throughput *= sample_transmittance;
diff --git a/crates/bevy_pbr/src/atmosphere/functions.wgsl b/crates/bevy_pbr/src/atmosphere/functions.wgsl
index daf66a6f8e60e..ebb8c4da984f7 100644
--- a/crates/bevy_pbr/src/atmosphere/functions.wgsl
+++ b/crates/bevy_pbr/src/atmosphere/functions.wgsl
@@ -43,6 +43,7 @@ const FRAC_3_16_PI: f32 = 0.0596831036594607509; // 3 / (16π)
 const FRAC_4_PI: f32 = 0.07957747154594767; // 1 / (4π)
 const ROOT_2: f32 = 1.41421356; // √2
 const EPSILON: f32 = 1.0; // 1 meter
+const MIN_EXTINCTION: vec3<f32> = vec3(1e-12);
 
 // During raymarching, each segment is sampled at a single point. This constant determines
 // where in the segment that sample is taken (0.0 = start, 0.5 = middle, 1.0 = end).
@@ -445,7 +446,7 @@ fn raymarch_atmosphere(
             sample_pos
         );
 
-        let s_int = (inscattering - inscattering * sample_transmittance) / extinction;
+        let s_int = (inscattering - inscattering * sample_transmittance) / max(extinction, MIN_EXTINCTION);
         result.inscattering += result.transmittance * s_int;
 
         result.transmittance *= sample_transmittance;
diff --git a/crates/bevy_pbr/src/atmosphere/multiscattering_lut.wgsl b/crates/bevy_pbr/src/atmosphere/multiscattering_lut.wgsl
index 03d55c9819961..8afe0c3be366f 100644
--- a/crates/bevy_pbr/src/atmosphere/multiscattering_lut.wgsl
+++ b/crates/bevy_pbr/src/atmosphere/multiscattering_lut.wgsl
@@ -7,7 +7,7 @@
             multiscattering_lut_uv_to_r_mu, sample_transmittance_lut,
             get_local_r, get_local_up, sample_density_lut, FRAC_4_PI,
             max_atmosphere_distance, rayleigh, henyey_greenstein,
-            zenith_azimuth_to_ray_dir,
+            zenith_azimuth_to_ray_dir, MIN_EXTINCTION
         },
         bruneton_functions::{
             distance_to_top_atmosphere_boundary, distance_to_bottom_atmosphere_boundary, ray_intersects_ground
@@ -109,7 +109,7 @@ fn sample_multiscattering_dir(r: f32, ray_dir: vec3<f32>, light_dir: vec3<f32>)
         optical_depth += sample_optical_depth;
 
         let ms = scattering;
-        let ms_int = (ms - ms * sample_transmittance) / extinction;
+        let ms_int = (ms - ms * sample_transmittance) / max(extinction, MIN_EXTINCTION);
         f_ms += throughput * ms_int;
 
         let mu_light = dot(light_dir, local_up);
@@ -117,7 +117,7 @@ fn sample_multiscattering_dir(r: f32, ray_dir: vec3<f32>, light_dir: vec3<f32>)
         let shadow_factor = transmittance_to_light * f32(!ray_intersects_ground(local_r, mu_light));
 
         let s = scattering * shadow_factor * FRAC_4_PI;
-        let s_int = (s - s * sample_transmittance) / extinction;
+        let s_int = (s - s * sample_transmittance) / max(extinction, MIN_EXTINCTION);
         l_2 += throughput * s_int;
 
         throughput *= sample_transmittance;