diff --git a/intern/cycles/blender/addon/properties.py b/intern/cycles/blender/addon/properties.py
index ca3936c10ad..329ed77424a 100644
--- a/intern/cycles/blender/addon/properties.py
+++ b/intern/cycles/blender/addon/properties.py
@@ -779,10 +779,32 @@ class CyclesRenderSettings(bpy.types.PropertyGroup):
         default=8,
     )
 
-    volume_unbiased: BoolProperty(
-        name="Unbiased",
-        description="If enabled, volume rendering converges to the correct result with sufficiently large numbers "
-        "of samples, but might appear noisier in the process",
+    volume_step_rate: FloatProperty(
+        name="Step Rate",
+        description="Globally adjust detail for volume rendering, on top of automatically estimated step size. "
+                    "Higher values reduce render time, lower values render with more detail",
+        default=1.0,
+        min=0.01, max=100.0, soft_min=0.1, soft_max=10.0, precision=2
+    )
+    volume_preview_step_rate: FloatProperty(
+        name="Step Rate",
+        description="Globally adjust detail for volume rendering, on top of automatically estimated step size. "
+                    "Higher values reduce render time, lower values render with more detail",
+        default=1.0,
+        min=0.01, max=100.0, soft_min=0.1, soft_max=10.0, precision=2
+    )
+    volume_max_steps: IntProperty(
+        name="Max Steps",
+        description="Maximum number of steps through the volume before giving up, "
+        "to avoid extremely long render times with big objects or small step sizes",
+        default=1024,
+        min=2, max=65536
+    )
+
+    volume_biased: BoolProperty(
+        name="Biased",
+        description="Default volume rendering uses null scattering, which is unbiased and has less artifacts, "
+        "but could be noisier. Biased option uses ray marching, with controls for steps size and max steps",
         default=False,
     )
 
@@ -1134,6 +1156,14 @@ class CyclesMaterialSettings(bpy.types.PropertyGroup):
         default='LINEAR',
     )
 
+    volume_step_rate: FloatProperty(
+        name="Step Rate",
+        description="Scale the distance between volume shader samples when rendering the volume "
+                    "(lower values give more accurate and detailed results, but also increased render time)",
+        default=1.0,
+        min=0.001, max=1000.0, soft_min=0.1, soft_max=10.0, precision=4
+    )
+
     @classmethod
     def register(cls):
         bpy.types.Material.cycles = PointerProperty(
@@ -1228,6 +1258,13 @@ class CyclesWorldSettings(bpy.types.PropertyGroup):
         items=enum_volume_interpolation,
         default='LINEAR',
     )
+    volume_step_size: FloatProperty(
+        name="Step Size",
+        description="Distance between volume shader samples when rendering the volume "
+                    "(lower values give more accurate and detailed results, but also increased render time)",
+        default=1.0,
+        min=0.0000001, max=100000.0, soft_min=0.1, soft_max=100.0, precision=4
+    )
 
     @classmethod
     def register(cls):
diff --git a/intern/cycles/blender/addon/ui.py b/intern/cycles/blender/addon/ui.py
index fc181744180..102bfb7cda3 100644
--- a/intern/cycles/blender/addon/ui.py
+++ b/intern/cycles/blender/addon/ui.py
@@ -584,8 +584,13 @@ class CYCLES_RENDER_PT_volumes(CyclesButtonsPanel, Panel):
         scene = context.scene
         cscene = scene.cycles
 
-        col = layout.column()
-        col.prop(cscene, "volume_unbiased", text="Unbiased")
+        col = layout.column(align=True)
+        col.prop(cscene, "volume_biased", text="Biased")
+        if cscene.volume_biased:
+            col.prop(cscene, "volume_step_rate", text="Step Rate Render")
+            col.prop(cscene, "volume_preview_step_rate", text="Viewport")
+
+            layout.prop(cscene, "volume_max_steps", text="Max Steps")
 
 
 class CYCLES_RENDER_PT_light_paths(CyclesButtonsPanel, Panel):
@@ -1836,6 +1841,8 @@ class CYCLES_WORLD_PT_settings_volume(CyclesButtonsPanel, Panel):
         sub = col.column()
         col.prop(cworld, "volume_sampling", text="Sampling")
         col.prop(cworld, "volume_interpolation", text="Interpolation")
+        if context.scene.cycles.volume_biased:
+            col.prop(cworld, "volume_step_size")
 
 
 class CYCLES_WORLD_PT_settings_light_group(CyclesButtonsPanel, Panel):
@@ -2008,6 +2015,8 @@ class CYCLES_MATERIAL_PT_settings_volume(CyclesButtonsPanel, Panel):
         sub = col.column()
         col.prop(cmat, "volume_sampling", text="Sampling")
         col.prop(cmat, "volume_interpolation", text="Interpolation")
+        if context.scene.cycles.volume_biased:
+            col.prop(cmat, "volume_step_rate")
 
     def draw(self, context):
         self.draw_shared(self, context, context.material)
diff --git a/intern/cycles/blender/shader.cpp b/intern/cycles/blender/shader.cpp
index 17634d334c3..efb2e7bcaa5 100644
--- a/intern/cycles/blender/shader.cpp
+++ b/intern/cycles/blender/shader.cpp
@@ -1625,6 +1625,7 @@ void BlenderSync::sync_materials(BL::Depsgraph &b_depsgraph, bool update_all)
       shader->set_use_bump_map_correction(get_boolean(cmat, "use_bump_map_correction"));
       shader->set_volume_sampling_method(get_volume_sampling(cmat));
       shader->set_volume_interpolation_method(get_volume_interpolation(cmat));
+      shader->set_volume_step_rate(get_float(cmat, "volume_step_rate"));
       shader->set_displacement_method(get_displacement_method(b_mat));
 
       shader->set_graph(std::move(graph));
@@ -1691,6 +1692,7 @@ void BlenderSync::sync_world(BL::Depsgraph &b_depsgraph, BL::SpaceView3D &b_v3d,
       PointerRNA cworld = RNA_pointer_get(&b_world.ptr, "cycles");
       shader->set_volume_sampling_method(get_volume_sampling(cworld));
       shader->set_volume_interpolation_method(get_volume_interpolation(cworld));
+      shader->set_volume_step_rate(get_float(cworld, "volume_step_size"));
     }
     else if (new_viewport_parameters.use_scene_world && b_world) {
       BackgroundNode *background = graph->create_node<BackgroundNode>();
diff --git a/intern/cycles/blender/sync.cpp b/intern/cycles/blender/sync.cpp
index 0e24821409e..239b234dea0 100644
--- a/intern/cycles/blender/sync.cpp
+++ b/intern/cycles/blender/sync.cpp
@@ -348,10 +348,15 @@ void BlenderSync::sync_integrator(BL::ViewLayer &b_view_layer,
   integrator->set_max_glossy_bounce(get_int(cscene, "glossy_bounces"));
   integrator->set_max_transmission_bounce(get_int(cscene, "transmission_bounces"));
   integrator->set_max_volume_bounce(get_int(cscene, "volume_bounces"));
-  integrator->set_volume_unbiased(get_boolean(cscene, "volume_unbiased"));
   integrator->set_transparent_min_bounce(get_int(cscene, "min_transparent_bounces"));
   integrator->set_transparent_max_bounce(get_int(cscene, "transparent_max_bounces"));
 
+  integrator->set_volume_ray_marching(get_boolean(cscene, "volume_biased"));
+  integrator->set_volume_max_steps(get_int(cscene, "volume_max_steps"));
+  const float volume_step_rate = (preview) ? get_float(cscene, "volume_preview_step_rate") :
+                                             get_float(cscene, "volume_step_rate");
+  integrator->set_volume_step_rate(volume_step_rate);
+
   integrator->set_caustics_reflective(get_boolean(cscene, "caustics_reflective"));
   integrator->set_caustics_refractive(get_boolean(cscene, "caustics_refractive"));
   integrator->set_filter_glossy(get_float(cscene, "blur_glossy"));
diff --git a/intern/cycles/blender/volume.cpp b/intern/cycles/blender/volume.cpp
index 1e19baa7c92..7b3e52fbd76 100644
--- a/intern/cycles/blender/volume.cpp
+++ b/intern/cycles/blender/volume.cpp
@@ -284,6 +284,7 @@ static void sync_volume_object(BL::BlendData &b_data,
 
   BL::VolumeRender b_render(b_volume.render());
 
+  volume->set_step_size(b_render.step_size());
   volume->set_object_space((b_render.space() == BL::VolumeRender::space_OBJECT));
 
   float velocity_scale = b_volume.velocity_scale();
diff --git a/intern/cycles/device/kernel.cpp b/intern/cycles/device/kernel.cpp
index bc8481789e6..ec59bbfd7ca 100644
--- a/intern/cycles/device/kernel.cpp
+++ b/intern/cycles/device/kernel.cpp
@@ -18,6 +18,7 @@ bool device_kernel_has_shading(DeviceKernel kernel)
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_RAYTRACE ||
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_MNEE ||
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME ||
+          kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING ||
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_SHADOW ||
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_DEDICATED_LIGHT ||
           kernel == DEVICE_KERNEL_SHADER_EVAL_DISPLACE ||
@@ -69,6 +70,8 @@ const char *device_kernel_as_string(DeviceKernel kernel)
       return "integrator_shade_surface_mnee";
     case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME:
       return "integrator_shade_volume";
+    case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING:
+      return "integrator_shade_volume_ray_marching";
     case DEVICE_KERNEL_INTEGRATOR_SHADE_DEDICATED_LIGHT:
       return "integrator_shade_dedicated_light";
     case DEVICE_KERNEL_INTEGRATOR_MEGAKERNEL:
diff --git a/intern/cycles/device/oneapi/device_impl.cpp b/intern/cycles/device/oneapi/device_impl.cpp
index 77536dae6dc..8363cb8bb2f 100644
--- a/intern/cycles/device/oneapi/device_impl.cpp
+++ b/intern/cycles/device/oneapi/device_impl.cpp
@@ -1284,6 +1284,7 @@ void OneapiDevice::get_adjusted_global_and_local_sizes(SyclQueue *queue,
     case DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_RAYTRACE:
     case DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_MNEE:
     case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME:
+    case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING:
     case DEVICE_KERNEL_INTEGRATOR_SHADE_SHADOW:
     case DEVICE_KERNEL_INTEGRATOR_SHADE_DEDICATED_LIGHT: {
       const bool device_is_simd8 =
diff --git a/intern/cycles/device/optix/device_impl.cpp b/intern/cycles/device/optix/device_impl.cpp
index 228761ea2fe..d28afbdd2c7 100644
--- a/intern/cycles/device/optix/device_impl.cpp
+++ b/intern/cycles/device/optix/device_impl.cpp
@@ -574,6 +574,10 @@ bool OptiXDevice::load_kernels(const uint kernel_features)
     group_descs[PG_RGEN_SHADE_VOLUME].raygen.module = optix_module;
     group_descs[PG_RGEN_SHADE_VOLUME].raygen.entryFunctionName =
         "__raygen__kernel_optix_integrator_shade_volume";
+    group_descs[PG_RGEN_SHADE_VOLUME_RAY_MARCHING].kind = OPTIX_PROGRAM_GROUP_KIND_RAYGEN;
+    group_descs[PG_RGEN_SHADE_VOLUME_RAY_MARCHING].raygen.module = optix_module;
+    group_descs[PG_RGEN_SHADE_VOLUME_RAY_MARCHING].raygen.entryFunctionName =
+        "__raygen__kernel_optix_integrator_shade_volume_ray_marching";
     group_descs[PG_RGEN_SHADE_SHADOW].kind = OPTIX_PROGRAM_GROUP_KIND_RAYGEN;
     group_descs[PG_RGEN_SHADE_SHADOW].raygen.module = optix_module;
     group_descs[PG_RGEN_SHADE_SHADOW].raygen.entryFunctionName =
diff --git a/intern/cycles/device/optix/device_impl.h b/intern/cycles/device/optix/device_impl.h
index f1db63feb9b..35626c1f940 100644
--- a/intern/cycles/device/optix/device_impl.h
+++ b/intern/cycles/device/optix/device_impl.h
@@ -32,6 +32,7 @@ enum {
   PG_RGEN_SHADE_SURFACE_RAYTRACE,
   PG_RGEN_SHADE_SURFACE_MNEE,
   PG_RGEN_SHADE_VOLUME,
+  PG_RGEN_SHADE_VOLUME_RAY_MARCHING,
   PG_RGEN_SHADE_SHADOW,
   PG_RGEN_SHADE_DEDICATED_LIGHT,
   PG_RGEN_EVAL_DISPLACE,
diff --git a/intern/cycles/device/optix/queue.cpp b/intern/cycles/device/optix/queue.cpp
index d86667e5230..be7cb6f51f8 100644
--- a/intern/cycles/device/optix/queue.cpp
+++ b/intern/cycles/device/optix/queue.cpp
@@ -124,6 +124,11 @@ bool OptiXDeviceQueue::enqueue(DeviceKernel kernel,
       pipeline = optix_device->pipelines[PIP_SHADE];
       sbt_params.raygenRecord = sbt_data_ptr + PG_RGEN_SHADE_VOLUME * sizeof(SbtRecord);
       break;
+    case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING:
+      pipeline = optix_device->pipelines[PIP_SHADE];
+      sbt_params.raygenRecord = sbt_data_ptr +
+                                PG_RGEN_SHADE_VOLUME_RAY_MARCHING * sizeof(SbtRecord);
+      break;
     case DEVICE_KERNEL_INTEGRATOR_SHADE_SHADOW:
       pipeline = optix_device->pipelines[PIP_SHADE];
       sbt_params.raygenRecord = sbt_data_ptr + PG_RGEN_SHADE_SHADOW * sizeof(SbtRecord);
diff --git a/intern/cycles/integrator/path_trace_work_gpu.cpp b/intern/cycles/integrator/path_trace_work_gpu.cpp
index af029e8ac42..9f5b147b52a 100644
--- a/intern/cycles/integrator/path_trace_work_gpu.cpp
+++ b/intern/cycles/integrator/path_trace_work_gpu.cpp
@@ -562,6 +562,7 @@ void PathTraceWorkGPU::enqueue_path_iteration(DeviceKernel kernel, const int num
     case DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_RAYTRACE:
     case DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_MNEE:
     case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME:
+    case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING:
     case DEVICE_KERNEL_INTEGRATOR_SHADE_DEDICATED_LIGHT: {
       /* Shading kernels with integrator state and render buffer. */
       const DeviceKernelArguments args(
@@ -570,7 +571,6 @@ void PathTraceWorkGPU::enqueue_path_iteration(DeviceKernel kernel, const int num
       queue_->enqueue(kernel, work_size, args);
       break;
     }
-
     default:
       LOG_FATAL << "Unhandled kernel " << device_kernel_as_string(kernel)
                 << " used for path iteration, should never happen.";
@@ -1272,6 +1272,7 @@ bool PathTraceWorkGPU::kernel_creates_shadow_paths(DeviceKernel kernel)
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_RAYTRACE ||
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_MNEE ||
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME ||
+          kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING ||
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_DEDICATED_LIGHT);
 }
 
diff --git a/intern/cycles/kernel/data_arrays.h b/intern/cycles/kernel/data_arrays.h
index 43c3d57e707..f4ad2abd8de 100644
--- a/intern/cycles/kernel/data_arrays.h
+++ b/intern/cycles/kernel/data_arrays.h
@@ -88,5 +88,6 @@ KERNEL_DATA_ARRAY(float, ies)
 KERNEL_DATA_ARRAY(KernelOctreeNode, volume_tree_nodes)
 KERNEL_DATA_ARRAY(KernelOctreeRoot, volume_tree_roots)
 KERNEL_DATA_ARRAY(int, volume_tree_root_ids)
+KERNEL_DATA_ARRAY(float, volume_step_size)
 
 #undef KERNEL_DATA_ARRAY
diff --git a/intern/cycles/kernel/data_template.h b/intern/cycles/kernel/data_template.h
index 74ad79f5db8..13eb1a6153d 100644
--- a/intern/cycles/kernel/data_template.h
+++ b/intern/cycles/kernel/data_template.h
@@ -208,7 +208,8 @@ KERNEL_STRUCT_MEMBER_DONT_SPECIALIZE
 KERNEL_STRUCT_MEMBER(integrator, int, blue_noise_sequence_length)
 /* Volume render. */
 KERNEL_STRUCT_MEMBER(integrator, int, use_volumes)
-KERNEL_STRUCT_MEMBER(integrator, int, volume_unbiased)
+KERNEL_STRUCT_MEMBER(integrator, int, volume_ray_marching)
+KERNEL_STRUCT_MEMBER(integrator, int, volume_max_steps)
 /* Shadow catcher. */
 KERNEL_STRUCT_MEMBER(integrator, int, has_shadow_catcher)
 /* Closure filter. */
@@ -231,7 +232,6 @@ KERNEL_STRUCT_MEMBER(integrator, int, use_guiding_mis_weights)
 /* Padding. */
 KERNEL_STRUCT_MEMBER(integrator, int, pad1)
 KERNEL_STRUCT_MEMBER(integrator, int, pad2)
-KERNEL_STRUCT_MEMBER(integrator, int, pad3)
 KERNEL_STRUCT_END(KernelIntegrator)
 
 /* SVM. For shader specialization. */
diff --git a/intern/cycles/kernel/device/gpu/kernel.h b/intern/cycles/kernel/device/gpu/kernel.h
index f237258925d..efcb06fb7f9 100644
--- a/intern/cycles/kernel/device/gpu/kernel.h
+++ b/intern/cycles/kernel/device/gpu/kernel.h
@@ -354,6 +354,21 @@ ccl_gpu_kernel(GPU_KERNEL_BLOCK_NUM_THREADS, GPU_KERNEL_MAX_REGISTERS)
 }
 ccl_gpu_kernel_postfix
 
+ccl_gpu_kernel(GPU_KERNEL_BLOCK_NUM_THREADS, GPU_KERNEL_MAX_REGISTERS)
+    ccl_gpu_kernel_signature(integrator_shade_volume_ray_marching,
+                             const ccl_global int *path_index_array,
+                             ccl_global float *render_buffer,
+                             const int work_size)
+{
+  const int global_index = ccl_gpu_global_id_x();
+
+  if (ccl_gpu_kernel_within_bounds(global_index, work_size)) {
+    const int state = (path_index_array) ? path_index_array[global_index] : global_index;
+    ccl_gpu_kernel_call(integrator_shade_volume_ray_marching(nullptr, state, render_buffer));
+  }
+}
+ccl_gpu_kernel_postfix
+
 ccl_gpu_kernel(GPU_KERNEL_BLOCK_NUM_THREADS, GPU_KERNEL_MAX_REGISTERS)
     ccl_gpu_kernel_signature(integrator_shade_dedicated_light,
                              const ccl_global int *path_index_array,
diff --git a/intern/cycles/kernel/device/oneapi/kernel.cpp b/intern/cycles/kernel/device/oneapi/kernel.cpp
index fbf2d0940c3..df477246a3a 100644
--- a/intern/cycles/kernel/device/oneapi/kernel.cpp
+++ b/intern/cycles/kernel/device/oneapi/kernel.cpp
@@ -460,6 +460,15 @@ bool oneapi_enqueue_kernel(KernelContext *kernel_context,
               kg, cgh, global_size, local_size, args, oneapi_kernel_integrator_shade_volume);
           break;
         }
+        case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING: {
+          oneapi_call(kg,
+                      cgh,
+                      global_size,
+                      local_size,
+                      args,
+                      oneapi_kernel_integrator_shade_volume_ray_marching);
+          break;
+        }
         case DEVICE_KERNEL_INTEGRATOR_SHADE_DEDICATED_LIGHT: {
           oneapi_call(kg,
                       cgh,
diff --git a/intern/cycles/kernel/device/optix/kernel_osl.cu b/intern/cycles/kernel/device/optix/kernel_osl.cu
index d32845ae3bc..faea22dda83 100644
--- a/intern/cycles/kernel/device/optix/kernel_osl.cu
+++ b/intern/cycles/kernel/device/optix/kernel_osl.cu
@@ -53,6 +53,15 @@ extern "C" __global__ void __raygen__kernel_optix_integrator_shade_volume()
   integrator_shade_volume(nullptr, path_index, kernel_params.render_buffer);
 }
 
+extern "C" __global__ void __raygen__kernel_optix_integrator_shade_volume_ray_marching()
+{
+  const int global_index = optixGetLaunchIndex().x;
+  const int path_index = (kernel_params.path_index_array) ?
+                             kernel_params.path_index_array[global_index] :
+                             global_index;
+  integrator_shade_volume_ray_marching(nullptr, path_index, kernel_params.render_buffer);
+}
+
 extern "C" __global__ void __raygen__kernel_optix_integrator_shade_shadow()
 {
   const int global_index = optixGetLaunchIndex().x;
diff --git a/intern/cycles/kernel/integrator/intersect_closest.h b/intern/cycles/kernel/integrator/intersect_closest.h
index 4977766d4bb..633abb4334c 100644
--- a/intern/cycles/kernel/integrator/intersect_closest.h
+++ b/intern/cycles/kernel/integrator/intersect_closest.h
@@ -233,7 +233,13 @@ ccl_device_forceinline void integrator_intersect_next_kernel(
     const int flags = (hit_surface) ? kernel_data_fetch(shaders, shader).flags : 0;
 
     if (!integrator_intersect_terminate(kg, state, flags)) {
-      integrator_path_next(state, current_kernel, DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME);
+      if (kernel_data.integrator.volume_ray_marching) {
+        integrator_path_next(
+            state, current_kernel, DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING);
+      }
+      else {
+        integrator_path_next(state, current_kernel, DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME);
+      }
     }
     else {
       integrator_path_terminate(kg, state, render_buffer, current_kernel);
diff --git a/intern/cycles/kernel/integrator/megakernel.h b/intern/cycles/kernel/integrator/megakernel.h
index 3e843cb0172..96b9791e134 100644
--- a/intern/cycles/kernel/integrator/megakernel.h
+++ b/intern/cycles/kernel/integrator/megakernel.h
@@ -76,6 +76,9 @@ ccl_device void integrator_megakernel(KernelGlobals kg,
         case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME:
           integrator_shade_volume(kg, state, render_buffer);
           break;
+        case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING:
+          integrator_shade_volume_ray_marching(kg, state, render_buffer);
+          break;
         case DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_RAYTRACE:
           integrator_shade_surface_raytrace(kg, state, render_buffer);
           break;
diff --git a/intern/cycles/kernel/integrator/shade_shadow.h b/intern/cycles/kernel/integrator/shade_shadow.h
index af5c0972f2d..daf2dbf78ed 100644
--- a/intern/cycles/kernel/integrator/shade_shadow.h
+++ b/intern/cycles/kernel/integrator/shade_shadow.h
@@ -94,7 +94,13 @@ ccl_device_inline void integrate_transparent_volume_shadow(KernelGlobals kg,
   /* `object` is only needed for light tree with light linking, it is irrelevant for shadow. */
   shader_setup_from_volume(shadow_sd, &ray, OBJECT_NONE);
 
-  volume_shadow_heterogeneous(kg, state, &ray, shadow_sd, throughput);
+  if (kernel_data.integrator.volume_ray_marching) {
+    const float step_size = volume_stack_step_size<true>(kg, state);
+    volume_shadow_ray_marching(kg, state, &ray, shadow_sd, throughput, step_size);
+  }
+  else {
+    volume_shadow_null_scattering(kg, state, &ray, shadow_sd, throughput);
+  }
 }
 #  endif
 
diff --git a/intern/cycles/kernel/integrator/shade_volume.h b/intern/cycles/kernel/integrator/shade_volume.h
index d4dc26ea72e..b12be528cf7 100644
--- a/intern/cycles/kernel/integrator/shade_volume.h
+++ b/intern/cycles/kernel/integrator/shade_volume.h
@@ -636,11 +636,11 @@ ccl_device Spectrum volume_transmittance(KernelGlobals kg,
 
 /* Compute the volumetric transmittance of the segment [ray->tmin, ray->tmax],
  * used for the shadow ray throughput. */
-ccl_device void volume_shadow_heterogeneous(KernelGlobals kg,
-                                            IntegratorShadowState state,
-                                            ccl_private Ray *ccl_restrict ray,
-                                            ccl_private ShaderData *ccl_restrict sd,
-                                            ccl_private Spectrum *ccl_restrict throughput)
+ccl_device void volume_shadow_null_scattering(KernelGlobals kg,
+                                              IntegratorShadowState state,
+                                              ccl_private Ray *ccl_restrict ray,
+                                              ccl_private ShaderData *ccl_restrict sd,
+                                              ccl_private Spectrum *ccl_restrict throughput)
 {
   /* Load random number state. */
   RNGState rng_state;
@@ -1891,14 +1891,14 @@ volume_direct_sample_method(KernelGlobals kg,
 }
 
 /* Shared function of integrating homogeneous and heterogeneous volume. */
-ccl_device void volume_integrate_shared(KernelGlobals kg,
-                                        const IntegratorState state,
-                                        const ccl_private Ray *ccl_restrict ray,
-                                        ccl_private ShaderData *ccl_restrict sd,
-                                        const ccl_private RNGState *rng_state,
-                                        ccl_global float *ccl_restrict render_buffer,
-                                        ccl_private LightSample *ls,
-                                        ccl_private VolumeIntegrateResult &result)
+ccl_device void volume_integrate_null_scattering(KernelGlobals kg,
+                                                 const IntegratorState state,
+                                                 const ccl_private Ray *ccl_restrict ray,
+                                                 ccl_private ShaderData *ccl_restrict sd,
+                                                 const ccl_private RNGState *rng_state,
+                                                 ccl_global float *ccl_restrict render_buffer,
+                                                 ccl_private LightSample *ls,
+                                                 ccl_private VolumeIntegrateResult &result)
 {
   PROFILING_INIT(kg, PROFILING_SHADE_VOLUME_INTEGRATE);
 
@@ -1944,6 +1944,426 @@ ccl_device void volume_integrate_shared(KernelGlobals kg,
   }
 }
 
+/* -------------------------------------------------------------------- */
+/** \name Ray Marching
+ * \{ */
+
+/* Determines the next shading position. */
+struct VolumeStep {
+  /* Shift starting point of all segments by a random amount to avoid banding artifacts due to
+   * biased ray marching with insufficient step size. */
+  float offset;
+
+  /* Step size taken at each marching step. */
+  float size;
+
+  /* Perform shading at this offset within a step, to integrate over the entire step segment. */
+  float shade_offset;
+
+  /* Maximal steps allowed between `ray->tmin` and `ray->tmax`. */
+  int max_steps;
+
+  /* Current active segment. */
+  Interval<float> t;
+};
+
+template<const bool shadow>
+ccl_device_forceinline void volume_step_init(KernelGlobals kg,
+                                             const ccl_private RNGState *rng_state,
+                                             const float object_step_size,
+                                             const float tmin,
+                                             const float tmax,
+                                             ccl_private VolumeStep *vstep)
+{
+  vstep->t.min = vstep->t.max = tmin;
+
+  if (object_step_size == FLT_MAX) {
+    /* Homogeneous volume. */
+    vstep->size = tmax - tmin;
+    vstep->shade_offset = 0.0f;
+    vstep->offset = 1.0f;
+    vstep->max_steps = 1;
+  }
+  else {
+    /* Heterogeneous volume. */
+    vstep->max_steps = kernel_data.integrator.volume_max_steps;
+    const float t = tmax - tmin;
+    float step_size = min(object_step_size, t);
+
+    if (t > vstep->max_steps * step_size) {
+      /* Increase step size to cover the whole ray segment. */
+      step_size = t / (float)vstep->max_steps;
+    }
+
+    vstep->size = step_size;
+    vstep->shade_offset = path_state_rng_1D(kg, rng_state, PRNG_VOLUME_SHADE_OFFSET);
+
+    if (shadow) {
+      /* For shadows we do not offset all segments, since the starting point is already a random
+       * distance inside the volume. It also appears to create banding artifacts for unknown
+       * reasons. */
+      vstep->offset = 1.0f;
+    }
+    else {
+      vstep->offset = path_state_rng_1D(kg, rng_state, PRNG_VOLUME_OFFSET);
+    }
+  }
+}
+
+ccl_device_inline bool volume_ray_marching_advance(const int step,
+                                                   const ccl_private Ray *ccl_restrict ray,
+                                                   ccl_private float3 *shade_P,
+                                                   ccl_private VolumeStep &vstep)
+{
+  if (vstep.t.max == ray->tmax) {
+    /* Reached the last segment. */
+    return false;
+  }
+
+  /* Advance to new position. */
+  vstep.t.min = vstep.t.max;
+  vstep.t.max = min(ray->tmax, ray->tmin + (step + vstep.offset) * vstep.size);
+  const float shade_t = mix(vstep.t.min, vstep.t.max, vstep.shade_offset);
+  *shade_P = ray->P + ray->D * shade_t;
+
+  return step < vstep.max_steps;
+}
+
+ccl_device void volume_shadow_ray_marching(KernelGlobals kg,
+                                           IntegratorShadowState state,
+                                           ccl_private Ray *ccl_restrict ray,
+                                           ccl_private ShaderData *ccl_restrict sd,
+                                           ccl_private Spectrum *ccl_restrict throughput,
+                                           const float object_step_size)
+{
+  /* Load random number state. */
+  RNGState rng_state;
+  shadow_path_state_rng_load(state, &rng_state);
+
+  /* For stochastic texture sampling. */
+  sd->lcg_state = lcg_state_init(
+      rng_state.rng_pixel, rng_state.rng_offset, rng_state.sample, 0xd9111870);
+
+  Spectrum tp = *throughput;
+
+  /* Prepare for stepping. */
+  VolumeStep vstep;
+  volume_step_init<true>(kg, &rng_state, object_step_size, ray->tmin, ray->tmax, &vstep);
+
+  /* compute extinction at the start */
+  Spectrum sum = zero_spectrum();
+  for (int step = 0; volume_ray_marching_advance(step, ray, &sd->P, vstep); step++) {
+    /* compute attenuation over segment */
+    const Spectrum sigma_t = volume_shader_eval_extinction<true>(kg, state, sd, PATH_RAY_SHADOW);
+    /* Compute `expf()` only for every Nth step, to save some calculations
+     * because `exp(a)*exp(b) = exp(a+b)`, also do a quick #VOLUME_THROUGHPUT_EPSILON
+     * check then. */
+    sum += (-sigma_t * vstep.t.length());
+    if ((step & 0x07) == 0) { /* TODO: Other interval? */
+      tp = *throughput * exp(sum);
+
+      /* stop if nearly all light is blocked */
+      if (reduce_max(tp) < VOLUME_THROUGHPUT_EPSILON) {
+        break;
+      }
+    }
+  }
+
+  if (vstep.t.max == ray->tmax) {
+    /* Update throughput in case we haven't done it above. */
+    tp = *throughput * exp(sum);
+  }
+
+  *throughput = tp;
+}
+
+struct VolumeRayMarchingState {
+  /* Random numbers for scattering. */
+  float rscatter;
+  float rchannel;
+
+  /* Multiple importance sampling. */
+  VolumeSampleMethod direct_sample_method;
+  bool use_mis;
+  float distance_pdf;
+  float equiangular_pdf;
+};
+
+ccl_device_inline void volume_ray_marching_state_init(
+    KernelGlobals kg,
+    const ccl_private RNGState *rng_state,
+    const VolumeSampleMethod direct_sample_method,
+    ccl_private VolumeRayMarchingState &vstate)
+{
+  vstate.rscatter = path_state_rng_1D(kg, rng_state, PRNG_VOLUME_SCATTER_DISTANCE);
+  vstate.rchannel = path_state_rng_1D(kg, rng_state, PRNG_VOLUME_COLOR_CHANNEL);
+
+  /* Multiple importance sampling: pick between equiangular and distance sampling strategy. */
+  vstate.direct_sample_method = direct_sample_method;
+  vstate.use_mis = (direct_sample_method == VOLUME_SAMPLE_MIS);
+  if (vstate.use_mis) {
+    if (vstate.rscatter < 0.5f) {
+      vstate.rscatter *= 2.0f;
+      vstate.direct_sample_method = VOLUME_SAMPLE_DISTANCE;
+    }
+    else {
+      vstate.rscatter = (vstate.rscatter - 0.5f) * 2.0f;
+      vstate.direct_sample_method = VOLUME_SAMPLE_EQUIANGULAR;
+    }
+  }
+  vstate.equiangular_pdf = 0.0f;
+  vstate.distance_pdf = 1.0f;
+}
+
+/* Returns true if we found the indirect scatter position within the current active ray segment. */
+ccl_device bool volume_sample_indirect_scatter_ray_marching(
+    const Spectrum transmittance,
+    const Spectrum channel_pdf,
+    const int channel,
+    const ccl_private ShaderData *ccl_restrict sd,
+    const ccl_private VolumeShaderCoefficients &ccl_restrict coeff,
+    const ccl_private Interval<float> &t,
+    ccl_private VolumeRayMarchingState &ccl_restrict vstate,
+    ccl_private VolumeIntegrateResult &ccl_restrict result)
+{
+  if (result.indirect_scatter) {
+    /* Already sampled indirect scatter position. */
+    return false;
+  }
+
+  /* If sampled distance does not go beyond the current segment, we have found the scatter
+   * position. Otherwise continue searching and accumulate the transmittance along the ray. */
+  const float sample_transmittance = volume_channel_get(transmittance, channel);
+  if (1.0f - vstate.rscatter >= sample_transmittance) {
+    /* Pick `sigma_t` from a random channel. */
+    const float sample_sigma_t = volume_channel_get(coeff.sigma_t, channel);
+
+    /* Generate the next distance using random walk, following exponential distribution
+     * p(dt) = sigma_t * exp(-sigma_t * dt). */
+    const float new_dt = -logf(1.0f - vstate.rscatter) / sample_sigma_t;
+    const float new_t = t.min + new_dt;
+
+    const Spectrum new_transmittance = volume_color_transmittance(coeff.sigma_t, new_dt);
+    /* PDF for density-based distance sampling is handled implicitly via
+     * transmittance / pdf = exp(-sigma_t * dt) / (sigma_t * exp(-sigma_t * dt)) = 1 / sigma_t. */
+    const float distance_pdf = dot(channel_pdf, coeff.sigma_t * new_transmittance);
+
+    if (vstate.distance_pdf * distance_pdf > VOLUME_SAMPLE_PDF_CUTOFF) {
+      /* Update throughput. */
+      result.indirect_scatter = true;
+      result.indirect_t = new_t;
+      result.indirect_throughput *= coeff.sigma_s * new_transmittance / distance_pdf;
+      if (vstate.direct_sample_method == VOLUME_SAMPLE_DISTANCE) {
+        vstate.distance_pdf *= distance_pdf;
+      }
+
+      volume_shader_copy_phases(&result.indirect_phases, sd);
+
+      return true;
+    }
+  }
+  else {
+    /* Update throughput. */
+    const float distance_pdf = dot(channel_pdf, transmittance);
+    result.indirect_throughput *= transmittance / distance_pdf;
+    if (vstate.direct_sample_method == VOLUME_SAMPLE_DISTANCE) {
+      vstate.distance_pdf *= distance_pdf;
+    }
+
+    /* Remap rscatter so we can reuse it and keep thing stratified. */
+    vstate.rscatter = 1.0f - (1.0f - vstate.rscatter) / sample_transmittance;
+  }
+
+  return false;
+}
+
+/* Find direct and indirect scatter positions. */
+ccl_device_forceinline void volume_ray_marching_step_scattering(
+    const ccl_private ShaderData *sd,
+    const ccl_private Ray *ray,
+    const ccl_private EquiangularCoefficients &equiangular_coeffs,
+    const ccl_private VolumeShaderCoefficients &ccl_restrict coeff,
+    const Spectrum transmittance,
+    const ccl_private Interval<float> &t,
+    ccl_private VolumeRayMarchingState &ccl_restrict vstate,
+    ccl_private VolumeIntegrateResult &ccl_restrict result)
+{
+  /* Pick random color channel for sampling the scatter distance. We use the Veach one-sample model
+   * with balance heuristic for the channels.
+   * Set `albedo` to 1 for the channel where extinction coefficient `sigma_t` is zero, to make sure
+   * that we sample a distance outside the current segment when that channel is picked, meaning
+   * light passes through without attenuation. */
+  const Spectrum albedo = safe_divide_color(coeff.sigma_s, coeff.sigma_t, 1.0f);
+  Spectrum channel_pdf;
+  const int channel = volume_sample_channel(
+      albedo, result.indirect_throughput, &vstate.rchannel, &channel_pdf);
+
+  /* Equiangular sampling for direct lighting. */
+  if (vstate.direct_sample_method == VOLUME_SAMPLE_EQUIANGULAR && !result.direct_scatter) {
+    if (t.contains(result.direct_t) && vstate.equiangular_pdf > VOLUME_SAMPLE_PDF_CUTOFF) {
+      const float new_dt = result.direct_t - t.min;
+      const Spectrum new_transmittance = volume_color_transmittance(coeff.sigma_t, new_dt);
+
+      result.direct_scatter = true;
+      result.direct_throughput *= coeff.sigma_s * new_transmittance / vstate.equiangular_pdf;
+      volume_shader_copy_phases(&result.direct_phases, sd);
+
+      /* Multiple importance sampling. */
+      if (vstate.use_mis) {
+        const float distance_pdf = vstate.distance_pdf *
+                                   dot(channel_pdf, coeff.sigma_t * new_transmittance);
+        const float mis_weight = 2.0f * power_heuristic(vstate.equiangular_pdf, distance_pdf);
+        result.direct_throughput *= mis_weight;
+      }
+    }
+    else {
+      result.direct_throughput *= transmittance;
+      vstate.distance_pdf *= dot(channel_pdf, transmittance);
+    }
+  }
+
+  /* Distance sampling for indirect and optional direct lighting. */
+  if (volume_sample_indirect_scatter_ray_marching(
+          transmittance, channel_pdf, channel, sd, coeff, t, vstate, result))
+  {
+    if (vstate.direct_sample_method == VOLUME_SAMPLE_DISTANCE) {
+      /* If using distance sampling for direct light, just copy parameters of indirect light
+       * since we scatter at the same point. */
+      result.direct_scatter = true;
+      result.direct_t = result.indirect_t;
+      result.direct_throughput = result.indirect_throughput;
+      volume_shader_copy_phases(&result.direct_phases, sd);
+
+      /* Multiple importance sampling. */
+      if (vstate.use_mis) {
+        const float equiangular_pdf = volume_equiangular_pdf(
+            ray, equiangular_coeffs, result.indirect_t);
+        const float mis_weight = power_heuristic(vstate.distance_pdf, equiangular_pdf);
+        result.direct_throughput *= 2.0f * mis_weight;
+      }
+    }
+  }
+}
+
+/* heterogeneous volume distance sampling: integrate stepping through the
+ * volume until we reach the end, get absorbed entirely, or run out of
+ * iterations. this does probabilistically scatter or get transmitted through
+ * for path tracing where we don't want to branch. */
+ccl_device_forceinline void volume_integrate_ray_marching(
+    KernelGlobals kg,
+    const IntegratorState state,
+    const ccl_private Ray *ccl_restrict ray,
+    ccl_private ShaderData *ccl_restrict sd,
+    const ccl_private RNGState *ccl_restrict rng_state,
+    ccl_global float *ccl_restrict render_buffer,
+    const float object_step_size,
+    ccl_private LightSample *ls,
+    ccl_private VolumeIntegrateResult &result)
+{
+  PROFILING_INIT(kg, PROFILING_SHADE_VOLUME_INTEGRATE);
+
+  EquiangularCoefficients equiangular_coeffs = {zero_float3(), {ray->tmin, ray->tmax}};
+  const VolumeSampleMethod direct_sample_method = volume_direct_sample_method(
+      kg, state, ray, sd, rng_state, &equiangular_coeffs, ls);
+
+  /* Prepare for stepping. */
+  VolumeStep vstep;
+  volume_step_init<false>(kg, rng_state, object_step_size, ray->tmin, ray->tmax, &vstep);
+
+  /* Initialize volume integration state. */
+  VolumeRayMarchingState vstate ccl_optional_struct_init;
+  volume_ray_marching_state_init(kg, rng_state, direct_sample_method, vstate);
+
+  /* Initialize volume integration result. */
+  const Spectrum throughput = INTEGRATOR_STATE(state, path, throughput);
+  result.direct_throughput = (vstate.direct_sample_method == VOLUME_SAMPLE_NONE) ?
+                                 zero_spectrum() :
+                                 throughput;
+  result.indirect_throughput = throughput;
+
+  /* Equiangular sampling: compute distance and PDF in advance. */
+  if (vstate.direct_sample_method == VOLUME_SAMPLE_EQUIANGULAR) {
+    result.direct_t = volume_equiangular_sample(
+        ray, equiangular_coeffs, vstate.rscatter, &vstate.equiangular_pdf);
+  }
+#  if defined(__PATH_GUIDING__)
+  result.direct_sample_method = vstate.direct_sample_method;
+#  endif
+
+#  ifdef __DENOISING_FEATURES__
+  const bool write_denoising_features = (INTEGRATOR_STATE(state, path, flag) &
+                                         PATH_RAY_DENOISING_FEATURES);
+  Spectrum accum_albedo = zero_spectrum();
+#  endif
+  Spectrum accum_emission = zero_spectrum();
+
+  for (int step = 0; volume_ray_marching_advance(step, ray, &sd->P, vstep); step++) {
+    /* compute segment */
+    VolumeShaderCoefficients coeff ccl_optional_struct_init;
+    if (volume_shader_sample(kg, state, sd, &coeff)) {
+      const int closure_flag = sd->flag;
+
+      /* Evaluate transmittance over segment. */
+      const float dt = vstep.t.length();
+      const Spectrum transmittance = (closure_flag & SD_EXTINCTION) ?
+                                         volume_color_transmittance(coeff.sigma_t, dt) :
+                                         one_spectrum();
+
+      /* Emission. */
+      if (closure_flag & SD_EMISSION) {
+        /* Only write emission before indirect light scatter position, since we terminate
+         * stepping at that point if we have already found a direct light scatter position. */
+        if (!result.indirect_scatter) {
+          const Spectrum emission = volume_emission_integrate(&coeff, closure_flag, dt);
+          accum_emission += result.indirect_throughput * emission;
+          guiding_record_volume_emission(kg, state, emission);
+        }
+      }
+
+      if (closure_flag & SD_SCATTER) {
+#  ifdef __DENOISING_FEATURES__
+        /* Accumulate albedo for denoising features. */
+        if (write_denoising_features && (closure_flag & SD_SCATTER)) {
+          const Spectrum albedo = safe_divide_color(coeff.sigma_s, coeff.sigma_t);
+          accum_albedo += result.indirect_throughput * albedo * (one_spectrum() - transmittance);
+        }
+#  endif
+
+        /* Scattering and absorption. */
+        volume_ray_marching_step_scattering(
+            sd, ray, equiangular_coeffs, coeff, transmittance, vstep.t, vstate, result);
+      }
+      else if (closure_flag & SD_EXTINCTION) {
+        /* Absorption only. */
+        result.indirect_throughput *= transmittance;
+        result.direct_throughput *= transmittance;
+      }
+
+      if (volume_integrate_should_stop(result)) {
+        break;
+      }
+    }
+  }
+
+  /* Write accumulated emission. */
+  if (!is_zero(accum_emission)) {
+    if (light_link_object_match(kg, light_link_receiver_forward(kg, state), sd->object)) {
+      film_write_volume_emission(
+          kg, state, accum_emission, render_buffer, object_lightgroup(kg, sd->object));
+    }
+  }
+
+#  ifdef __DENOISING_FEATURES__
+  /* Write denoising features. */
+  if (write_denoising_features) {
+    film_write_denoising_features_volume(
+        kg, state, accum_albedo, result.indirect_scatter, render_buffer);
+  }
+#  endif /* __DENOISING_FEATURES__ */
+}
+
+/** \} */
+
 /* Path tracing: sample point on light and evaluate light shader, then
  * queue shadow ray to be traced. */
 ccl_device_forceinline void integrate_volume_direct_light(
@@ -2199,38 +2619,15 @@ ccl_device_forceinline bool integrate_volume_phase_scatter(
   return true;
 }
 
-/* get the volume attenuation and emission over line segment defined by
- * ray, with the assumption that there are no surfaces blocking light
- * between the endpoints. distance sampling is used to decide if we will
- * scatter or not. */
-ccl_device VolumeIntegrateEvent volume_integrate(KernelGlobals kg,
-                                                 IntegratorState state,
-                                                 ccl_private Ray *ccl_restrict ray,
-                                                 ccl_global float *ccl_restrict render_buffer)
+ccl_device_inline VolumeIntegrateEvent
+volume_integrate_event(KernelGlobals kg,
+                       IntegratorState state,
+                       const ccl_private Ray *ccl_restrict ray,
+                       ccl_private ShaderData *sd,
+                       const ccl_private RNGState *rng_state,
+                       ccl_private LightSample &ls,
+                       ccl_private VolumeIntegrateResult &result)
 {
-  if (integrator_state_volume_stack_is_empty(kg, state)) {
-    return VOLUME_PATH_ATTENUATED;
-  }
-
-  ShaderData sd;
-  /* FIXME: `object` is used for light linking. We read the bottom of the stack for simplicity, but
-   * this does not work for overlapping volumes. */
-  shader_setup_from_volume(&sd, ray, INTEGRATOR_STATE_ARRAY(state, volume_stack, 0, object));
-
-  /* Load random number state. */
-  RNGState rng_state;
-  path_state_rng_load(state, &rng_state);
-
-  /* For stochastic texture sampling. */
-  sd.lcg_state = lcg_state_init(
-      rng_state.rng_pixel, rng_state.rng_offset, rng_state.sample, 0x15b4f88d);
-
-  LightSample ls ccl_optional_struct_init;
-
-  /* TODO: expensive to zero closures? */
-  VolumeIntegrateResult result = {};
-  volume_integrate_shared(kg, state, ray, &sd, &rng_state, render_buffer, &ls, result);
-
 #  if defined(__PATH_GUIDING__) && PATH_GUIDING_LEVEL >= 1
   /* The current path throughput which is used later to calculate per-segment throughput. */
   const float3 initial_throughput = INTEGRATOR_STATE(state, path, throughput);
@@ -2267,10 +2664,10 @@ ccl_device VolumeIntegrateEvent volume_integrate(KernelGlobals kg,
         const float3 transmittance_weight = spectrum_to_rgb(
             safe_divide_color(result.indirect_throughput, initial_throughput));
         guiding_record_volume_transmission(kg, state, transmittance_weight);
-        guiding_record_volume_segment(kg, state, direct_P, sd.wi);
+        guiding_record_volume_segment(kg, state, direct_P, sd->wi);
         guiding_generated_new_segment = true;
         unlit_throughput = result.indirect_throughput / continuation_probability;
-        rand_phase_guiding = path_state_rng_1D(kg, &rng_state, PRNG_VOLUME_PHASE_GUIDING_DISTANCE);
+        rand_phase_guiding = path_state_rng_1D(kg, rng_state, PRNG_VOLUME_PHASE_GUIDING_DISTANCE);
       }
       else if (result.direct_sample_method == VOLUME_SAMPLE_EQUIANGULAR) {
         /* If the direct scatter event is generated using VOLUME_SAMPLE_EQUIANGULAR the direct
@@ -2287,7 +2684,7 @@ ccl_device VolumeIntegrateEvent volume_integrate(KernelGlobals kg,
         unlit_throughput /= scatterEval;
         unlit_throughput *= continuation_probability;
         rand_phase_guiding = path_state_rng_1D(
-            kg, &rng_state, PRNG_VOLUME_PHASE_GUIDING_EQUIANGULAR);
+            kg, rng_state, PRNG_VOLUME_PHASE_GUIDING_EQUIANGULAR);
       }
 #    endif
 #    if PATH_GUIDING_LEVEL >= 4
@@ -2302,8 +2699,8 @@ ccl_device VolumeIntegrateEvent volume_integrate(KernelGlobals kg,
     result.direct_throughput /= continuation_probability;
     integrate_volume_direct_light(kg,
                                   state,
-                                  &sd,
-                                  &rng_state,
+                                  sd,
+                                  rng_state,
                                   direct_P,
                                   &result.direct_phases,
 #  if defined(__PATH_GUIDING__)
@@ -2330,28 +2727,28 @@ ccl_device VolumeIntegrateEvent volume_integrate(KernelGlobals kg,
   INTEGRATOR_STATE_WRITE(state, path, throughput) = result.indirect_throughput;
 
   if (result.indirect_scatter) {
-    sd.P = ray->P + result.indirect_t * ray->D;
+    sd->P = ray->P + result.indirect_t * ray->D;
 
 #  if defined(__PATH_GUIDING__)
     if ((kernel_data.kernel_features & KERNEL_FEATURE_PATH_GUIDING)) {
 #    if PATH_GUIDING_LEVEL >= 1
       if (!guiding_generated_new_segment) {
-        guiding_record_volume_segment(kg, state, sd.P, sd.wi);
+        guiding_record_volume_segment(kg, state, sd->P, sd->wi);
       }
 #    endif
 #    if PATH_GUIDING_LEVEL >= 4
       /* If the direct scatter event was generated using VOLUME_SAMPLE_EQUIANGULAR we need to
        * initialize the guiding distribution at the indirect scatter position. */
       if (result.direct_sample_method == VOLUME_SAMPLE_EQUIANGULAR) {
-        rand_phase_guiding = path_state_rng_1D(kg, &rng_state, PRNG_VOLUME_PHASE_GUIDING_DISTANCE);
+        rand_phase_guiding = path_state_rng_1D(kg, rng_state, PRNG_VOLUME_PHASE_GUIDING_DISTANCE);
         volume_shader_prepare_guiding(
-            kg, state, rand_phase_guiding, sd.P, ray->D, &result.indirect_phases);
+            kg, state, rand_phase_guiding, sd->P, ray->D, &result.indirect_phases);
       }
 #    endif
     }
 #  endif
 
-    if (integrate_volume_phase_scatter(kg, state, &sd, ray, &rng_state, &result.indirect_phases)) {
+    if (integrate_volume_phase_scatter(kg, state, sd, ray, rng_state, &result.indirect_phases)) {
       return VOLUME_PATH_SCATTERED;
     }
     return VOLUME_PATH_MISSED;
@@ -2365,27 +2762,99 @@ ccl_device VolumeIntegrateEvent volume_integrate(KernelGlobals kg,
   return VOLUME_PATH_ATTENUATED;
 }
 
+/* get the volume attenuation and emission over line segment defined by
+ * ray, with the assumption that there are no surfaces blocking light
+ * between the endpoints. distance sampling is used to decide if we will
+ * scatter or not. */
+ccl_device VolumeIntegrateEvent volume_integrate(KernelGlobals kg,
+                                                 IntegratorState state,
+                                                 ccl_private Ray *ccl_restrict ray,
+                                                 ccl_global float *ccl_restrict render_buffer)
+{
+  kernel_assert(!kernel_data.integrator.volume_ray_marching);
+
+  if (integrator_state_volume_stack_is_empty(kg, state)) {
+    return VOLUME_PATH_ATTENUATED;
+  }
+
+  ShaderData sd;
+  /* FIXME: `object` is used for light linking. We read the bottom of the stack for simplicity, but
+   * this does not work for overlapping volumes. */
+  shader_setup_from_volume(&sd, ray, INTEGRATOR_STATE_ARRAY(state, volume_stack, 0, object));
+
+  /* Load random number state. */
+  RNGState rng_state;
+  path_state_rng_load(state, &rng_state);
+
+  /* For stochastic texture sampling. */
+  sd.lcg_state = lcg_state_init(
+      rng_state.rng_pixel, rng_state.rng_offset, rng_state.sample, 0x15b4f88d);
+
+  LightSample ls ccl_optional_struct_init;
+
+  /* TODO: expensive to zero closures? */
+  VolumeIntegrateResult result = {};
+  volume_integrate_null_scattering(kg, state, ray, &sd, &rng_state, render_buffer, &ls, result);
+
+  return volume_integrate_event(kg, state, ray, &sd, &rng_state, ls, result);
+}
+
+ccl_device VolumeIntegrateEvent
+volume_integrate_ray_marching(KernelGlobals kg,
+                              IntegratorState state,
+                              ccl_private Ray *ccl_restrict ray,
+                              ccl_global float *ccl_restrict render_buffer)
+{
+  kernel_assert(kernel_data.integrator.volume_ray_marching);
+
+  if (integrator_state_volume_stack_is_empty(kg, state)) {
+    return VOLUME_PATH_ATTENUATED;
+  }
+
+  ShaderData sd;
+  /* FIXME: `object` is used for light linking. We read the bottom of the stack for simplicity, but
+   * this does not work for overlapping volumes. */
+  shader_setup_from_volume(&sd, ray, INTEGRATOR_STATE_ARRAY(state, volume_stack, 0, object));
+
+  /* Load random number state. */
+  RNGState rng_state;
+  path_state_rng_load(state, &rng_state);
+
+  /* For stochastic texture sampling. */
+  sd.lcg_state = lcg_state_init(
+      rng_state.rng_pixel, rng_state.rng_offset, rng_state.sample, 0x15b4f88d);
+
+  LightSample ls ccl_optional_struct_init;
+
+  /* TODO: expensive to zero closures? */
+  VolumeIntegrateResult result = {};
+
+  const float step_size = volume_stack_step_size<false>(kg, state);
+  volume_integrate_ray_marching(
+      kg, state, ray, &sd, &rng_state, render_buffer, step_size, &ls, result);
+
+  return volume_integrate_event(kg, state, ray, &sd, &rng_state, ls, result);
+}
+
 #endif
 
-ccl_device void integrator_shade_volume(KernelGlobals kg,
-                                        IntegratorState state,
-                                        ccl_global float *ccl_restrict render_buffer)
+#ifdef __VOLUME__
+ccl_device_inline void integrator_shade_volume_setup(KernelGlobals kg,
+                                                     IntegratorState state,
+                                                     ccl_private Ray *ccl_restrict ray,
+                                                     ccl_private Intersection *ccl_restrict isect)
 {
   PROFILING_INIT(kg, PROFILING_SHADE_VOLUME_SETUP);
 
-#ifdef __VOLUME__
   /* Setup shader data. */
-  Ray ray ccl_optional_struct_init;
-  integrator_state_read_ray(state, &ray);
-
-  Intersection isect ccl_optional_struct_init;
-  integrator_state_read_isect(state, &isect);
+  integrator_state_read_ray(state, ray);
+  integrator_state_read_isect(state, isect);
 
   /* Set ray length to current segment. */
-  ray.tmax = (isect.prim != PRIM_NONE) ? isect.t : FLT_MAX;
+  ray->tmax = (isect->prim != PRIM_NONE) ? isect->t : FLT_MAX;
 
   /* Clean volume stack for background rays. */
-  if (isect.prim == PRIM_NONE) {
+  if (isect->prim == PRIM_NONE) {
     volume_stack_clean(kg, state);
   }
 
@@ -2393,35 +2862,70 @@ ccl_device void integrator_shade_volume(KernelGlobals kg,
   if (INTEGRATOR_STATE(state, path, bounce) == 0) {
     INTEGRATOR_STATE_WRITE(state, path, flag) |= PATH_RAY_VOLUME_PRIMARY_TRANSMIT;
   }
+}
+#endif
 
-  const VolumeIntegrateEvent event = volume_integrate(kg, state, &ray, render_buffer);
+template<DeviceKernel volume_kernel>
+ccl_device_inline void integrator_next_kernel_after_shade_volume(
+    KernelGlobals kg,
+    const IntegratorState state,
+    ccl_global float *ccl_restrict render_buffer,
+    const ccl_private Intersection *ccl_restrict isect,
+    const VolumeIntegrateEvent event)
+{
   if (event == VOLUME_PATH_MISSED) {
     /* End path. */
-    integrator_path_terminate(kg, state, render_buffer, DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME);
+    integrator_path_terminate(kg, state, render_buffer, volume_kernel);
     return;
   }
 
   if (event == VOLUME_PATH_ATTENUATED) {
     /* Continue to background, light or surface. */
-    integrator_intersect_next_kernel_after_volume<DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME>(
-        kg, state, &isect, render_buffer);
+    integrator_intersect_next_kernel_after_volume<volume_kernel>(kg, state, isect, render_buffer);
     return;
   }
 
-#  ifdef __SHADOW_LINKING__
-  if (shadow_linking_schedule_intersection_kernel<DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME>(kg,
-                                                                                         state))
-  {
+#ifdef __SHADOW_LINKING__
+  if (shadow_linking_schedule_intersection_kernel<volume_kernel>(kg, state)) {
     return;
   }
-#  endif /* __SHADOW_LINKING__ */
+#endif /* __SHADOW_LINKING__ */
 
   kernel_assert(event == VOLUME_PATH_SCATTERED);
 
   /* Queue intersect_closest kernel. */
-  integrator_path_next(
-      state, DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME, DEVICE_KERNEL_INTEGRATOR_INTERSECT_CLOSEST);
+  integrator_path_next(state, volume_kernel, DEVICE_KERNEL_INTEGRATOR_INTERSECT_CLOSEST);
+}
+
+ccl_device void integrator_shade_volume(KernelGlobals kg,
+                                        IntegratorState state,
+                                        ccl_global float *ccl_restrict render_buffer)
+{
+#ifdef __VOLUME__
+  Ray ray ccl_optional_struct_init;
+  Intersection isect ccl_optional_struct_init;
+  integrator_shade_volume_setup(kg, state, &ray, &isect);
+
+  const VolumeIntegrateEvent event = volume_integrate(kg, state, &ray, render_buffer);
+  integrator_next_kernel_after_shade_volume<DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME>(
+      kg, state, render_buffer, &isect, event);
+
 #endif /* __VOLUME__ */
 }
 
+ccl_device void integrator_shade_volume_ray_marching(KernelGlobals kg,
+                                                     IntegratorState state,
+                                                     ccl_global float *ccl_restrict render_buffer)
+{
+#ifdef __VOLUME__
+  Ray ray ccl_optional_struct_init;
+  Intersection isect ccl_optional_struct_init;
+  integrator_shade_volume_setup(kg, state, &ray, &isect);
+
+  const VolumeIntegrateEvent event = volume_integrate_ray_marching(kg, state, &ray, render_buffer);
+  integrator_next_kernel_after_shade_volume<DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING>(
+      kg, state, render_buffer, &isect, event);
+
+#endif /* __VOLUME__ */
+}
 CCL_NAMESPACE_END
diff --git a/intern/cycles/kernel/integrator/volume_stack.h b/intern/cycles/kernel/integrator/volume_stack.h
index f45b6dc37da..35b0cc9f462 100644
--- a/intern/cycles/kernel/integrator/volume_stack.h
+++ b/intern/cycles/kernel/integrator/volume_stack.h
@@ -185,6 +185,28 @@ ccl_device_inline bool volume_is_homogeneous(KernelGlobals kg, const IntegratorG
   return false;
 }
 
+template<const bool shadow, typename IntegratorGenericState>
+ccl_device float volume_stack_step_size(KernelGlobals kg, const IntegratorGenericState state)
+{
+  kernel_assert(kernel_data.integrator.volume_ray_marching);
+
+  float step_size = FLT_MAX;
+
+  for (int i = 0;; i++) {
+    const VolumeStack entry = volume_stack_read<shadow>(state, i);
+    if (entry.shader == SHADER_NONE) {
+      break;
+    }
+
+    if (!volume_is_homogeneous(kg, entry)) {
+      const float object_step_size = kernel_data_fetch(volume_step_size, entry.object);
+      step_size = fminf(object_step_size, step_size);
+    }
+  }
+
+  return step_size;
+}
+
 enum VolumeSampleMethod {
   VOLUME_SAMPLE_NONE = 0,
   VOLUME_SAMPLE_DISTANCE = (1 << 0),
diff --git a/intern/cycles/kernel/types.h b/intern/cycles/kernel/types.h
index 55de43cedb5..629eeaa97b0 100644
--- a/intern/cycles/kernel/types.h
+++ b/intern/cycles/kernel/types.h
@@ -297,6 +297,8 @@ enum PathTraceDimension {
   PRNG_VOLUME_SHADE_OFFSET = 7,
   PRNG_VOLUME_PHASE_GUIDING_DISTANCE = 8,
   PRNG_VOLUME_PHASE_GUIDING_EQUIANGULAR = 9,
+  PRNG_VOLUME_COLOR_CHANNEL = 4,
+  PRNG_VOLUME_OFFSET = 6,
 
   /* Subsurface random walk bounces */
   PRNG_SUBSURFACE_BSDF = 0,
@@ -1853,6 +1855,7 @@ enum DeviceKernel : int {
   DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_RAYTRACE,
   DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_MNEE,
   DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME,
+  DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING,
   DEVICE_KERNEL_INTEGRATOR_SHADE_SHADOW,
   DEVICE_KERNEL_INTEGRATOR_SHADE_DEDICATED_LIGHT,
   DEVICE_KERNEL_INTEGRATOR_MEGAKERNEL,
diff --git a/intern/cycles/scene/devicescene.cpp b/intern/cycles/scene/devicescene.cpp
index 27aa20ba1d1..d3630014d15 100644
--- a/intern/cycles/scene/devicescene.cpp
+++ b/intern/cycles/scene/devicescene.cpp
@@ -56,7 +56,8 @@ DeviceScene::DeviceScene(Device *device)
       ies_lights(device, "ies", MEM_GLOBAL),
       volume_tree_nodes(device, "volume_tree_nodes", MEM_GLOBAL),
       volume_tree_roots(device, "volume_tree_roots", MEM_GLOBAL),
-      volume_tree_root_ids(device, "volume_tree_root_ids", MEM_GLOBAL)
+      volume_tree_root_ids(device, "volume_tree_root_ids", MEM_GLOBAL),
+      volume_step_size(device, "volume_step_size", MEM_GLOBAL)
 {
   memset((void *)&data, 0, sizeof(data));
 }
diff --git a/intern/cycles/scene/devicescene.h b/intern/cycles/scene/devicescene.h
index 9aaa7bbc22e..9dbd2ecba62 100644
--- a/intern/cycles/scene/devicescene.h
+++ b/intern/cycles/scene/devicescene.h
@@ -89,6 +89,7 @@ class DeviceScene {
   device_vector<KernelOctreeNode> volume_tree_nodes;
   device_vector<KernelOctreeRoot> volume_tree_roots;
   device_vector<int> volume_tree_root_ids;
+  device_vector<float> volume_step_size;
 
   KernelData data;
 
diff --git a/intern/cycles/scene/integrator.cpp b/intern/cycles/scene/integrator.cpp
index 2610b4dafe2..eef15d0757a 100644
--- a/intern/cycles/scene/integrator.cpp
+++ b/intern/cycles/scene/integrator.cpp
@@ -57,7 +57,9 @@ NODE_DEFINE(Integrator)
   SOCKET_FLOAT(ao_distance, "AO Distance", FLT_MAX);
   SOCKET_FLOAT(ao_additive_factor, "AO Additive Factor", 0.0f);
 
-  SOCKET_BOOLEAN(volume_unbiased, "Unbiased", false);
+  SOCKET_BOOLEAN(volume_ray_marching, "Biased", false);
+  SOCKET_INT(volume_max_steps, "Volume Max Steps", 1024);
+  SOCKET_FLOAT(volume_step_rate, "Volume Step Rate", 1.0f);
 
   static NodeEnum guiding_distribution_enum;
   guiding_distribution_enum.insert("PARALLAX_AWARE_VMM", GUIDING_TYPE_PARALLAX_AWARE_VMM);
@@ -226,7 +228,8 @@ void Integrator::device_update(Device *device, DeviceScene *dscene, Scene *scene
     }
   }
 
-  kintegrator->volume_unbiased = volume_unbiased;
+  kintegrator->volume_ray_marching = volume_ray_marching;
+  kintegrator->volume_max_steps = volume_max_steps;
 
   kintegrator->caustics_reflective = caustics_reflective;
   kintegrator->caustics_refractive = caustics_refractive;
diff --git a/intern/cycles/scene/integrator.h b/intern/cycles/scene/integrator.h
index bb2ed42c792..0672250a53b 100644
--- a/intern/cycles/scene/integrator.h
+++ b/intern/cycles/scene/integrator.h
@@ -41,7 +41,9 @@ class Integrator : public Node {
   NODE_SOCKET_API(float, ao_distance)
   NODE_SOCKET_API(float, ao_additive_factor)
 
-  NODE_SOCKET_API(bool, volume_unbiased)
+  NODE_SOCKET_API(bool, volume_ray_marching)
+  NODE_SOCKET_API(int, volume_max_steps)
+  NODE_SOCKET_API(float, volume_step_rate)
 
   NODE_SOCKET_API(bool, use_guiding);
   NODE_SOCKET_API(bool, deterministic_guiding);
diff --git a/intern/cycles/scene/object.cpp b/intern/cycles/scene/object.cpp
index f96920a6220..5112270fb07 100644
--- a/intern/cycles/scene/object.cpp
+++ b/intern/cycles/scene/object.cpp
@@ -289,6 +289,97 @@ uint Object::visibility_for_tracing() const
   return SHADOW_CATCHER_OBJECT_VISIBILITY(is_shadow_catcher, visibility & PATH_RAY_ALL_VISIBILITY);
 }
 
+float Object::compute_volume_step_size() const
+{
+  if (geometry->is_light()) {
+    /* World volume. */
+    assert(static_cast<const Light *>(geometry)->get_light_type() == LIGHT_BACKGROUND);
+    for (const Node *node : geometry->get_used_shaders()) {
+      const Shader *shader = static_cast<const Shader *>(node);
+      if (shader->has_volume) {
+        return shader->get_volume_step_rate();
+      }
+    }
+    assert(false);
+    return FLT_MAX;
+  }
+
+  if (!geometry->is_mesh() && !geometry->is_volume()) {
+    return FLT_MAX;
+  }
+
+  Mesh *mesh = static_cast<Mesh *>(geometry);
+
+  if (!mesh->has_volume) {
+    return FLT_MAX;
+  }
+
+  /* Compute step rate from shaders. */
+  float step_rate = FLT_MAX;
+
+  for (Node *node : mesh->get_used_shaders()) {
+    Shader *shader = static_cast<Shader *>(node);
+    if (shader->has_volume) {
+      if (shader->has_volume_spatial_varying || shader->has_volume_attribute_dependency) {
+        step_rate = fminf(shader->get_volume_step_rate(), step_rate);
+      }
+    }
+  }
+
+  if (step_rate == FLT_MAX) {
+    return FLT_MAX;
+  }
+
+  /* Compute step size from voxel grids. */
+  float step_size = FLT_MAX;
+
+  if (geometry->is_volume()) {
+    Volume *volume = static_cast<Volume *>(geometry);
+
+    for (Attribute &attr : volume->attributes.attributes) {
+      if (attr.element == ATTR_ELEMENT_VOXEL) {
+        ImageHandle &handle = attr.data_voxel();
+        const ImageMetaData &metadata = handle.metadata();
+        if (metadata.byte_size == 0) {
+          continue;
+        }
+
+        /* User specified step size. */
+        float voxel_step_size = volume->get_step_size();
+
+        if (voxel_step_size == 0.0f) {
+          /* Auto detect step size.
+           * Step size is transformed from voxel to world space. */
+          Transform voxel_tfm = tfm;
+          if (metadata.use_transform_3d) {
+            voxel_tfm = tfm * transform_inverse(metadata.transform_3d);
+          }
+          voxel_step_size = reduce_min(fabs(transform_direction(&voxel_tfm, one_float3())));
+        }
+        else if (volume->get_object_space()) {
+          /* User specified step size in object space. */
+          const float3 size = make_float3(voxel_step_size, voxel_step_size, voxel_step_size);
+          voxel_step_size = reduce_min(fabs(transform_direction(&tfm, size)));
+        }
+
+        if (voxel_step_size > 0.0f) {
+          step_size = fminf(voxel_step_size, step_size);
+        }
+      }
+    }
+  }
+
+  if (step_size == FLT_MAX) {
+    /* Fall back to 1/10th of bounds for procedural volumes. */
+    assert(bounds.valid());
+    step_size = 0.1f * average(bounds.size());
+  }
+
+  step_size *= step_rate;
+
+  return step_size;
+}
+
 int Object::get_device_index() const
 {
   return index;
@@ -689,6 +780,7 @@ void ObjectManager::device_update(Device *device,
     dscene->object_motion_pass.tag_realloc();
     dscene->object_motion.tag_realloc();
     dscene->object_flag.tag_realloc();
+    dscene->volume_step_size.tag_realloc();
 
     /* If objects are added to the scene or deleted, the object indices might change, so we need to
      * update the root indices of the volume octrees. */
@@ -730,6 +822,7 @@ void ObjectManager::device_update(Device *device,
         dscene->object_motion_pass.tag_modified();
         dscene->object_motion.tag_modified();
         dscene->object_flag.tag_modified();
+        dscene->volume_step_size.tag_modified();
       }
 
       /* Update world object index. */
@@ -800,6 +893,9 @@ void ObjectManager::device_update_flags(Device * /*unused*/,
   bool has_volume_objects = false;
   for (Object *object : scene->objects) {
     if (object->geometry->has_volume) {
+      /* If the bounds are not valid it is not always possible to calculate the volume step, and
+       * the step size is not needed for the displacement. So, delay calculation of the volume
+       * step size until the final bounds are known. */
       if (bounds_valid) {
         volume_objects.push_back(object);
       }
@@ -852,7 +948,6 @@ void ObjectManager::device_update_flags(Device * /*unused*/,
 
   /* Copy object flag. */
   dscene->object_flag.copy_to_device();
-
   dscene->object_flag.clear_modified();
 }
 
diff --git a/intern/cycles/scene/object.h b/intern/cycles/scene/object.h
index 170fcef1e6b..026bc8dbcc2 100644
--- a/intern/cycles/scene/object.h
+++ b/intern/cycles/scene/object.h
@@ -112,6 +112,9 @@ class Object : public Node {
   /* Returns the index that is used in the kernel for this object. */
   int get_device_index() const;
 
+  /* Compute step size from attributes, shaders, transforms. */
+  float compute_volume_step_size() const;
+
   /* Check whether this object can be used as light-emissive. */
   bool usable_as_light() const;
 
diff --git a/intern/cycles/scene/shader.cpp b/intern/cycles/scene/shader.cpp
index 69c6b1f3dd3..5fc1b0d6f3f 100644
--- a/intern/cycles/scene/shader.cpp
+++ b/intern/cycles/scene/shader.cpp
@@ -72,6 +72,8 @@ NODE_DEFINE(Shader)
               volume_interpolation_method_enum,
               VOLUME_INTERPOLATION_LINEAR);
 
+  SOCKET_FLOAT(volume_step_rate, "Volume Step Rate", 1.0f);
+
   static NodeEnum displacement_method_enum;
   displacement_method_enum.insert("bump", DISPLACE_BUMP);
   displacement_method_enum.insert("true", DISPLACE_TRUE);
@@ -101,6 +103,7 @@ Shader::Shader() : Node(get_node_type())
   has_volume_spatial_varying = false;
   has_volume_attribute_dependency = false;
   has_volume_connected = false;
+  prev_volume_step_rate = 0.0f;
   has_light_path_node = false;
 
   emission_estimate = zero_float3();
@@ -390,9 +393,10 @@ void Shader::tag_update(Scene *scene)
     scene->procedural_manager->tag_update();
   }
 
-  if (has_volume != prev_has_volume) {
+  if (has_volume != prev_has_volume || volume_step_rate != prev_volume_step_rate) {
     scene->geometry_manager->need_flags_update = true;
     scene->object_manager->need_flags_update = true;
+    prev_volume_step_rate = volume_step_rate;
   }
 
   if (has_volume || prev_has_volume) {
diff --git a/intern/cycles/scene/shader.h b/intern/cycles/scene/shader.h
index 7b323c30119..f14ddaa4c00 100644
--- a/intern/cycles/scene/shader.h
+++ b/intern/cycles/scene/shader.h
@@ -83,10 +83,13 @@ class Shader : public Node {
   NODE_SOCKET_API(bool, use_bump_map_correction)
   NODE_SOCKET_API(VolumeSampling, volume_sampling_method)
   NODE_SOCKET_API(int, volume_interpolation_method)
+  NODE_SOCKET_API(float, volume_step_rate)
 
   /* displacement */
   NODE_SOCKET_API(DisplacementMethod, displacement_method)
 
+  float prev_volume_step_rate;
+
   /* synchronization */
   bool need_update_uvs;
   bool need_update_attribute;
diff --git a/intern/cycles/scene/volume.cpp b/intern/cycles/scene/volume.cpp
index 8aa38bf1710..f7b602c235b 100644
--- a/intern/cycles/scene/volume.cpp
+++ b/intern/cycles/scene/volume.cpp
@@ -6,6 +6,7 @@
 #include "scene/attribute.h"
 #include "scene/background.h"
 #include "scene/image_vdb.h"
+#include "scene/integrator.h"
 #include "scene/light.h"
 #include "scene/object.h"
 #include "scene/scene.h"
@@ -34,6 +35,7 @@ NODE_DEFINE(Volume)
 {
   NodeType *type = NodeType::add("volume", create, NodeType::NONE, Mesh::get_node_type());
 
+  SOCKET_FLOAT(step_size, "Step Size", 0.0f);
   SOCKET_BOOLEAN(object_space, "Object Space", false);
   SOCKET_FLOAT(velocity_scale, "Velocity Scale", 1.0f);
 
@@ -42,6 +44,7 @@ NODE_DEFINE(Volume)
 
 Volume::Volume() : Mesh(get_node_type(), Geometry::VOLUME)
 {
+  step_size = 0.0f;
   object_space = false;
 }
 
@@ -156,9 +159,12 @@ class VolumeMeshBuilder {
 
   void create_mesh(vector<float3> &vertices,
                    vector<int> &indices,
-                   const float face_overlap_avoidance);
+                   const float face_overlap_avoidance,
+                   const bool ray_marching);
 
-  void generate_vertices_and_quads(vector<int3> &vertices_is, vector<QuadData> &quads);
+  void generate_vertices_and_quads(vector<int3> &vertices_is,
+                                   vector<QuadData> &quads,
+                                   const bool ray_marching);
 
   void convert_object_space(const vector<int3> &vertices,
                             vector<float3> &out_vertices,
@@ -207,23 +213,96 @@ void VolumeMeshBuilder::add_padding(const int pad_size)
 
 void VolumeMeshBuilder::create_mesh(vector<float3> &vertices,
                                     vector<int> &indices,
-                                    const float face_overlap_avoidance)
+                                    const float face_overlap_avoidance,
+                                    const bool ray_marching)
 {
   /* We create vertices in index space (is), and only convert them to object
    * space when done. */
   vector<int3> vertices_is;
   vector<QuadData> quads;
 
-  generate_vertices_and_quads(vertices_is, quads);
+  generate_vertices_and_quads(vertices_is, quads, ray_marching);
 
   convert_object_space(vertices_is, vertices, face_overlap_avoidance);
 
   convert_quads_to_tris(quads, indices);
 }
 
-void VolumeMeshBuilder::generate_vertices_and_quads(vector<ccl::int3> &vertices_is,
-                                                    vector<QuadData> &quads)
+static bool is_non_empty_leaf(const openvdb::MaskGrid::TreeType &tree, const openvdb::Coord coord)
 {
+  const auto *leaf_node = tree.probeLeaf(coord);
+  return (leaf_node && !leaf_node->isEmpty());
+}
+
+void VolumeMeshBuilder::generate_vertices_and_quads(vector<ccl::int3> &vertices_is,
+                                                    vector<QuadData> &quads,
+                                                    const bool ray_marching)
+{
+  if (ray_marching) {
+    /* Make sure we only have leaf nodes in the tree, as tiles are not handled by this algorithm */
+    topology_grid->tree().voxelizeActiveTiles();
+
+    const openvdb::MaskGrid::TreeType &tree = topology_grid->tree();
+    tree.evalLeafBoundingBox(bbox);
+
+    const int3 resolution = make_int3(bbox.dim().x(), bbox.dim().y(), bbox.dim().z());
+
+    unordered_map<size_t, int> used_verts;
+    for (auto iter = tree.cbeginLeaf(); iter; ++iter) {
+      if (iter->isEmpty()) {
+        continue;
+      }
+      openvdb::CoordBBox leaf_bbox = iter->getNodeBoundingBox();
+      /* +1 to convert from exclusive to include bounds. */
+      leaf_bbox.max() = leaf_bbox.max().offsetBy(1);
+      int3 min = make_int3(leaf_bbox.min().x(), leaf_bbox.min().y(), leaf_bbox.min().z());
+      int3 max = make_int3(leaf_bbox.max().x(), leaf_bbox.max().y(), leaf_bbox.max().z());
+      int3 corners[8] = {
+          make_int3(min[0], min[1], min[2]),
+          make_int3(max[0], min[1], min[2]),
+          make_int3(max[0], max[1], min[2]),
+          make_int3(min[0], max[1], min[2]),
+          make_int3(min[0], min[1], max[2]),
+          make_int3(max[0], min[1], max[2]),
+          make_int3(max[0], max[1], max[2]),
+          make_int3(min[0], max[1], max[2]),
+      };
+      /* Only create a quad if on the border between an active and an inactive leaf.
+       *
+       * We verify that a leaf exists by probing a coordinate that is at its center,
+       * to do so we compute the center of the current leaf and offset this coordinate
+       * by the size of a leaf in each direction.
+       */
+      static const int LEAF_DIM = openvdb::MaskGrid::TreeType::LeafNodeType::DIM;
+      auto center = leaf_bbox.min() + openvdb::Coord(LEAF_DIM / 2);
+      if (!is_non_empty_leaf(tree, openvdb::Coord(center.x() - LEAF_DIM, center.y(), center.z())))
+      {
+        create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_X_MIN);
+      }
+      if (!is_non_empty_leaf(tree, openvdb::Coord(center.x() + LEAF_DIM, center.y(), center.z())))
+      {
+        create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_X_MAX);
+      }
+      if (!is_non_empty_leaf(tree, openvdb::Coord(center.x(), center.y() - LEAF_DIM, center.z())))
+      {
+        create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_Y_MIN);
+      }
+      if (!is_non_empty_leaf(tree, openvdb::Coord(center.x(), center.y() + LEAF_DIM, center.z())))
+      {
+        create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_Y_MAX);
+      }
+      if (!is_non_empty_leaf(tree, openvdb::Coord(center.x(), center.y(), center.z() - LEAF_DIM)))
+      {
+        create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_Z_MIN);
+      }
+      if (!is_non_empty_leaf(tree, openvdb::Coord(center.x(), center.y(), center.z() + LEAF_DIM)))
+      {
+        create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_Z_MAX);
+      }
+    }
+    return;
+  }
+
   bbox = topology_grid->evalActiveVoxelBoundingBox();
 
   const int3 resolution = make_int3(bbox.dim().x(), bbox.dim().y(), bbox.dim().z());
@@ -524,7 +603,8 @@ void GeometryManager::create_volume_mesh(const Scene *scene, Volume *volume, Pro
   /* Create mesh. */
   vector<float3> vertices;
   vector<int> indices;
-  builder.create_mesh(vertices, indices, face_overlap_avoidance);
+  const bool ray_marching = scene->integrator->get_volume_ray_marching();
+  builder.create_mesh(vertices, indices, face_overlap_avoidance, ray_marching);
 
   volume->reserve_mesh(vertices.size(), indices.size() / 3);
   volume->used_shaders.clear();
@@ -570,6 +650,11 @@ void VolumeManager::tag_update()
 /* Remove changed object from the list of octrees and tag for rebuild. */
 void VolumeManager::tag_update(const Object *object, uint32_t flag)
 {
+  if (object_octrees_.empty()) {
+    /* Volume object is not in the octree, can happen when using ray marching. */
+    return;
+  }
+
   if (flag & ObjectManager::VISIBILITY_MODIFIED) {
     tag_update();
   }
@@ -1013,12 +1098,52 @@ std::string VolumeManager::visualize_octree(const char *filename) const
   return filename_full;
 }
 
+void VolumeManager::update_step_size(const Scene *scene, DeviceScene *dscene) const
+{
+  assert(scene->integrator->get_volume_ray_marching());
+
+  if (!dscene->volume_step_size.is_modified() && last_algorithm == RAY_MARCHING) {
+    return;
+  }
+
+  if (dscene->volume_step_size.need_realloc()) {
+    dscene->volume_step_size.alloc(scene->objects.size());
+  }
+
+  float *volume_step_size = dscene->volume_step_size.data();
+
+  for (const Object *object : scene->objects) {
+    const Geometry *geom = object->get_geometry();
+    if (!geom->has_volume) {
+      continue;
+    }
+
+    volume_step_size[object->index] = scene->integrator->get_volume_step_rate() *
+                                      object->compute_volume_step_size();
+  }
+
+  dscene->volume_step_size.copy_to_device();
+  dscene->volume_step_size.clear_modified();
+}
+
 void VolumeManager::device_update(Device *device,
                                   DeviceScene *dscene,
                                   const Scene *scene,
                                   Progress &progress)
 {
-  if (need_rebuild_) {
+  if (scene->integrator->get_volume_ray_marching()) {
+    /* No need to update octree for ray marching. */
+    if (last_algorithm == NULL_SCATTERING) {
+      dscene->volume_tree_nodes.free();
+      dscene->volume_tree_roots.free();
+      dscene->volume_tree_root_ids.free();
+    }
+    update_step_size(scene, dscene);
+    last_algorithm = RAY_MARCHING;
+    return;
+  }
+
+  if (need_rebuild_ || last_algorithm == RAY_MARCHING) {
     /* Data needed for volume shader evaluation. */
     device->const_copy_to("data", &dscene->data, sizeof(dscene->data));
 
@@ -1039,6 +1164,11 @@ void VolumeManager::device_update(Device *device,
     LOG_DEBUG << "Octree visualization has been written to " << visualize_octree("octree.py");
     update_visualization_ = false;
   }
+
+  if (last_algorithm == RAY_MARCHING) {
+    dscene->volume_step_size.free();
+  }
+  last_algorithm = NULL_SCATTERING;
 }
 
 void VolumeManager::device_free(DeviceScene *dscene)
@@ -1046,6 +1176,7 @@ void VolumeManager::device_free(DeviceScene *dscene)
   dscene->volume_tree_nodes.free();
   dscene->volume_tree_roots.free();
   dscene->volume_tree_root_ids.free();
+  dscene->volume_step_size.free();
 }
 
 VolumeManager::~VolumeManager()
diff --git a/intern/cycles/scene/volume.h b/intern/cycles/scene/volume.h
index 46fdfbfe7fb..984ffa7d416 100644
--- a/intern/cycles/scene/volume.h
+++ b/intern/cycles/scene/volume.h
@@ -17,12 +17,19 @@ CCL_NAMESPACE_BEGIN
 class Object;
 class Octree;
 
+enum VolumeRenderingAlgorithm {
+  NULL_SCATTERING,
+  RAY_MARCHING,
+  NONE,
+};
+
 class Volume : public Mesh {
  public:
   NODE_DECLARE
 
   Volume();
 
+  NODE_SOCKET_API(float, step_size)
   NODE_SOCKET_API(bool, object_space)
   NODE_SOCKET_API(float, velocity_scale)
 
@@ -72,6 +79,9 @@ class VolumeManager {
    * `filename`, which is a Python script that can be run inside Blender. */
   std::string visualize_octree(const char *filename) const;
 
+  /* Step size for ray marching. */
+  void update_step_size(const Scene *, DeviceScene *) const;
+
   /* One octree per object per shader. */
   std::map<std::pair<const Object *, const Shader *>, std::shared_ptr<Octree>> object_octrees_;
 
@@ -81,6 +91,8 @@ class VolumeManager {
   int num_octree_nodes_;
   int num_octree_roots_;
 
+  VolumeRenderingAlgorithm last_algorithm = NONE;
+
 #ifdef WITH_OPENVDB
   /* Create SDF grid for mesh volumes, to determine whether a certain point is in the
    * interior of the mesh. This reduces evaluation time needed for heterogeneous volume. */
diff --git a/scripts/startup/bl_ui/properties_data_volume.py b/scripts/startup/bl_ui/properties_data_volume.py
index d6ddf6958e1..00f4af7792a 100644
--- a/scripts/startup/bl_ui/properties_data_volume.py
+++ b/scripts/startup/bl_ui/properties_data_volume.py
@@ -128,6 +128,8 @@ class DATA_PT_volume_render(DataButtonsPanel, Panel):
         col.prop(render, "space")
 
         if scene.render.engine == 'CYCLES':
+            col.prop(render, "step_size")
+
             col = layout.column(align=True)
             col.prop(render, "clipping")
 
diff --git a/source/blender/makesrna/intern/rna_volume.cc b/source/blender/makesrna/intern/rna_volume.cc
index 33ea7c28e00..f9708564f59 100644
--- a/source/blender/makesrna/intern/rna_volume.cc
+++ b/source/blender/makesrna/intern/rna_volume.cc
@@ -541,6 +541,17 @@ static void rna_def_volume_render(BlenderRNA *brna)
       prop, "Space", "Specify volume density and step size in object or world space");
   RNA_def_property_update(prop, 0, "rna_Volume_update_display");
 
+  prop = RNA_def_property(srna, "step_size", PROP_FLOAT, PROP_DISTANCE);
+  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
+  RNA_def_property_range(prop, 0.0, FLT_MAX);
+  RNA_def_property_ui_range(prop, 0.0, 100.0, 1, 3);
+  RNA_def_property_ui_text(prop,
+                           "Step Size",
+                           "Distance between volume samples. Lower values render more detail at "
+                           "the cost of performance. If set to zero, the step size is "
+                           "automatically determined based on voxel size.");
+  RNA_def_property_update(prop, 0, "rna_Volume_update_display");
+
   prop = RNA_def_property(srna, "clipping", PROP_FLOAT, PROP_NONE);
   RNA_def_property_float_sdna(prop, nullptr, "clipping");
   RNA_def_property_range(prop, 0.0, 1.0);
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_camera_and_volume.png b/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_camera_and_volume.png
new file mode 100644
index 00000000000..ea7d1369900
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_camera_and_volume.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7572b66f9ca7f20a5ac1a4697d58eb73b7d26c1df0ed8dc7232de64c17f8642e
+size 31066
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_camera_only.png b/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_camera_only.png
new file mode 100644
index 00000000000..ce95ba9dd26
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_camera_only.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:adbfb1fa2d3aa39ab56977cefa648c96c2ee1161c5975c7786296fb26c70634d
+size 32198
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_nested.png b/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_nested.png
new file mode 100644
index 00000000000..6f63af7cca3
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_nested.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:b19317985ecbeab0101313e38d897db1a0a7f5163b21ff97a564277cc12672f3
+size 34080
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_simple.png b/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_simple.png
new file mode 100644
index 00000000000..3e2fd788c08
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_simple.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:dcc9ae25fb61bccc21a2cd18eda8011bff44f59bf886d01f0b16ff44f38d8225
+size 32407
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/implicit_volume.png b/tests/files/render/volume/cycles_ray_marching_renders/implicit_volume.png
new file mode 100644
index 00000000000..dd615aefbdf
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/implicit_volume.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:97eae519b8285147f3f9a799886e4973d9d432a50b45497502d023c714b8f7de
+size 45924
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/overlapping_different_anisotropy.png b/tests/files/render/volume/cycles_ray_marching_renders/overlapping_different_anisotropy.png
new file mode 100644
index 00000000000..d7db939d43b
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/overlapping_different_anisotropy.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:02441326ba412e6fa6b103cd18dd7429203fbdc70758ccde585844070c2ba4d1
+size 37744
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/overlapping_octrees.png b/tests/files/render/volume/cycles_ray_marching_renders/overlapping_octrees.png
new file mode 100644
index 00000000000..5a43614a408
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/overlapping_octrees.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c3308fcdf3b240a2fad27d08a7164fffc3c6e1a2e9616fb401af8ce406d7d15d
+size 25479
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/principled_absorption.png b/tests/files/render/volume/cycles_ray_marching_renders/principled_absorption.png
new file mode 100644
index 00000000000..6f3bf47367c
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/principled_absorption.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:b0a51cbd9bd51945acd917b137c67f5a04eb4c49ebc34b2a9321ff673833563d
+size 29659
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/principled_bsdf_interior.png b/tests/files/render/volume/cycles_ray_marching_renders/principled_bsdf_interior.png
new file mode 100644
index 00000000000..3cce76d91a5
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/principled_bsdf_interior.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a786c1532d647a03fb809fd18649bca09ee6e3537a1792ef76672723d215247d
+size 29335
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_absobtion.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_absobtion.png
new file mode 100644
index 00000000000..189a252ede2
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_absobtion.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:1731b6e31e4e0210c38921dabb96e800563d24d209dcaebab73396fb59ec22be
+size 23589
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_deep_stack.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_deep_stack.png
new file mode 100644
index 00000000000..605b33f9836
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_deep_stack.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a0a0d4666a11ff6af51a1204220d2b0fdf7c13c0c4bb4c438ed57cc3af695070
+size 4141
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_edge_fireflies.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_edge_fireflies.png
new file mode 100644
index 00000000000..29f6858a82b
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_edge_fireflies.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:dcfb133dbc7e7240c2e1fd787186e5935d4cd7e93c00d7e4a027519ea5f5f713
+size 3298
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_instance.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_instance.png
new file mode 100644
index 00000000000..d03d6ddd2c9
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_instance.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7d993a35b298d5a21616387c84edb0c194c354e7cd58f0a71be44af8329c4ebe
+size 28354
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_light_path.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_light_path.png
new file mode 100644
index 00000000000..52587087999
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_light_path.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ab7c1a2289edb13d746b740a6c001ca482da7d38fb3a7d4b65c9abb269ecaeaf
+size 23088
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_output_absorption.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_output_absorption.png
new file mode 100644
index 00000000000..0e5c87f4869
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_output_absorption.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:cc15cddd5027cf590328a6d21a9f9dd4cfb52ead19fc2a7052e6928a90024fd7
+size 23903
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_output_mix_volume.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_output_mix_volume.png
new file mode 100644
index 00000000000..a38ab36b5b5
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_output_mix_volume.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a528037e96290f5306330d1fe656b11859339818e999ebc51b979feef2fbb9c9
+size 25603
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_output_scatter.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_output_scatter.png
new file mode 100644
index 00000000000..5ef5959c128
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_output_scatter.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:de2cf1e9d414b4eb18c6895c774356d411b33c857935b3a91e32769d9395bda6
+size 26116
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_output_surface.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_output_surface.png
new file mode 100644
index 00000000000..6d076157aa2
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_output_surface.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f3d2801e019ca3a2136124b5c62b6486227514e7e78a34e4d58a20a4f7b3de31
+size 25025
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_overlap.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_overlap.png
new file mode 100644
index 00000000000..b212ee7aefa
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_overlap.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:65ffd7efcabedd956c5dfb51eb3abb148804c3731181274148920e46e9d043aa
+size 25998
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter.png
new file mode 100644
index 00000000000..6c4adf9b9a7
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0ff643c9c84624544a2d05b2ed070adbab245ce9f00dc83adf718389c4cfce51
+size 26109
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_albedo.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_albedo.png
new file mode 100644
index 00000000000..d4c976f5694
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_albedo.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:11d937b2f5724e88dcac3555c63dfc8e28c51fc5545dc862d971a7c6c8814708
+size 32745
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_draine.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_draine.png
new file mode 100644
index 00000000000..dbfc1a4431e
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_draine.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:592e2d7c80b0ebc4db8985fcd22d18ab4123b841a799ba7d95679e4277e60fb8
+size 26681
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_fournier-forand.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_fournier-forand.png
new file mode 100644
index 00000000000..d57b83e832b
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_fournier-forand.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:2dc3e38e17e7edb7719de4d2a5c5808a3eaceaebe98cb5a64cd96688960d9549
+size 36306
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_mie.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_mie.png
new file mode 100644
index 00000000000..93d21eca68c
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_mie.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:e59f9a1176a5d278406ee1cab45ec04f791163743de859a85091a0c50e75ac7f
+size 33926
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_mie_small_particles.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_mie_small_particles.png
new file mode 100644
index 00000000000..500ef4b5576
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_mie_small_particles.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:78810cc672185f94c1cd83586d605a8ebb05bd52f96cb0043470cac469799a1d
+size 34232
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_rayleigh.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_rayleigh.png
new file mode 100644
index 00000000000..f332ced7af7
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_rayleigh.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fe4b1b20c4c812af686614898e90a8a40e9270152a8e86e14590e71e5f8bcefc
+size 43458
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_step_offset.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_step_offset.png
new file mode 100644
index 00000000000..28b2059fc5e
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_step_offset.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c9046550e6fb8ebcaa36d4ee4e5e8648a2a4ce36edda564e321652a086e0bde3
+size 9410
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_transparent_shadow.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_transparent_shadow.png
new file mode 100644
index 00000000000..36085dfe5c0
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_transparent_shadow.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:3b49e7828249be2f02ba7033999919ac22d01534cdab1b6adcabf9f9832c0009
+size 32327
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_two_shaders.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_two_shaders.png
new file mode 100644
index 00000000000..29a43e25b08
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_two_shaders.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:86c1954c00369c5177a1b033bbfb71c1f601eef4b8e5bbd2ccf9f50a3bc76209
+size 27522
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_zero_extinction_channel.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_zero_extinction_channel.png
new file mode 100644
index 00000000000..f16cce2f009
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_zero_extinction_channel.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:88247095729b93c07dac899809633d8333e5b8bf52523321e75047e90cbd8ad4
+size 42210
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/world_volume.png b/tests/files/render/volume/cycles_ray_marching_renders/world_volume.png
new file mode 100644
index 00000000000..040c9fb6b8c
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/world_volume.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9ebb133a43bf16fc85199ee45d7dfc702a85e5b7d7f5255536340779ff9a8a6f
+size 31393
diff --git a/tests/python/cycles_render_tests.py b/tests/python/cycles_render_tests.py
index 1f8b69007c0..d5b8ef72d80 100644
--- a/tests/python/cycles_render_tests.py
+++ b/tests/python/cycles_render_tests.py
@@ -128,29 +128,39 @@ BLOCKLIST_GPU = [
 
 
 class CyclesReport(render_report.Report):
-    def __init__(self, title, output_dir, oiiotool, device=None, blocklist=[], osl=False):
+    def __init__(self, title, output_dir, oiiotool, device=None, blocklist=[], osl=False, ray_marching=False):
         # Split device name in format "<device_type>[-<RT>]" into individual
         # tokens, setting the RT suffix to an empty string if its not specified.
         self.device, suffix = (device.split("-") + [""])[:2]
         self.use_hwrt = (suffix == "RT")
         self.osl = osl
+        self.ray_marching = ray_marching
 
         variation = self.device
         if suffix:
             variation += ' ' + suffix
         if self.osl:
             variation += ' OSL'
+        if ray_marching:
+            variation += ' Ray Marching'
 
         super().__init__(title, output_dir, oiiotool, variation, blocklist)
 
+        self.set_pixelated(True)
+        self.set_reference_dir("cycles_renders")
+        if device == 'CPU':
+            self.set_compare_engine('eevee')
+        else:
+            self.set_compare_engine('cycles', 'CPU')
+
     def _get_render_arguments(self, arguments_cb, filepath, base_output_filepath):
-        return arguments_cb(filepath, base_output_filepath, self.use_hwrt, self.osl)
+        return arguments_cb(filepath, base_output_filepath, self.use_hwrt, self.osl, self.ray_marching)
 
     def _get_arguments_suffix(self):
         return ['--', '--cycles-device', self.device] if self.device else []
 
 
-def get_arguments(filepath, output_filepath, use_hwrt=False, osl=False):
+def get_arguments(filepath, output_filepath, use_hwrt=False, osl=False, ray_marching=False):
     dirname = os.path.dirname(filepath)
     basedir = os.path.dirname(dirname)
     subject = os.path.basename(dirname)
@@ -191,6 +201,9 @@ def get_arguments(filepath, output_filepath, use_hwrt=False, osl=False):
     if osl:
         args.extend(["--python-expr", "import bpy; bpy.context.scene.cycles.shading_system = True"])
 
+    if ray_marching:
+        args.extend(["--python-expr", "import bpy; bpy.context.scene.cycles.volume_biased = True"])
+
     if subject == 'bake':
         args.extend(['--python', os.path.join(basedir, "util", "render_bake.py")])
     elif subject == 'denoise_animation':
@@ -215,6 +228,13 @@ def create_argparse():
     return parser
 
 
+def test_volume_ray_marching(args, device, blocklist):
+    # Default volume rendering algorithm is null scattering, but we also want to test ray marching
+    report = CyclesReport('Cycles', args.outdir, args.oiiotool, device, blocklist, args.osl == 'all', ray_marching=True)
+    report.set_reference_dir("cycles_ray_marching_renders")
+    return report.run(args.testdir, args.blender, get_arguments, batch=args.batch)
+
+
 def main():
     parser = create_argparse()
     args = parser.parse_args()
@@ -243,12 +263,6 @@ def main():
         blocklist += BLOCKLIST_METAL
 
     report = CyclesReport('Cycles', args.outdir, args.oiiotool, device, blocklist, args.osl == 'all')
-    report.set_pixelated(True)
-    report.set_reference_dir("cycles_renders")
-    if device == 'CPU':
-        report.set_compare_engine('eevee')
-    else:
-        report.set_compare_engine('cycles', 'CPU')
 
     # Increase threshold for motion blur, see #78777.
     #
@@ -286,6 +300,9 @@ def main():
 
     ok = report.run(args.testdir, args.blender, get_arguments, batch=args.batch)
 
+    if (test_dir_name == 'volume'):
+        ok = ok and test_volume_ray_marching(args, device, blocklist)
+
     sys.exit(not ok)