diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl
index 8feb80f30..d21a68497 100644
--- a/examples/fluid/pipe_flow_2d.jl
+++ b/examples/fluid/pipe_flow_2d.jl
@@ -152,11 +152,12 @@ ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=100)
 saving_callback = SolutionSavingCallback(dt=0.02, prefix="")
+particle_shifting = ParticleShiftingCallback()
 
 extra_callback = nothing
 
 callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(),
-                        ParticleShiftingCallback(), extra_callback)
+                        particle_shifting, extra_callback)
 
 sol = solve(ode, RDPK3SpFSAL35(),
             abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl
index 6739f80f9..30e589e1c 100644
--- a/src/schemes/boundary/open_boundary/boundary_zones.jl
+++ b/src/schemes/boundary/open_boundary/boundary_zones.jl
@@ -7,7 +7,8 @@ struct OutFlow end
 @doc raw"""
     BoundaryZone(; plane, plane_normal, density, particle_spacing,
                  initial_condition=nothing, extrude_geometry=nothing,
-                 open_boundary_layers::Integer, boundary_type=BidirectionalFlow())
+                 open_boundary_layers::Integer, boundary_type=BidirectionalFlow(),
+                 average_inflow_velocity=true)
 
 Boundary zone for [`OpenBoundarySPHSystem`](@ref).
 
@@ -54,6 +55,14 @@ There are three ways to specify the actual shape of the boundary zone:
 - `extrude_geometry=nothing`: 1D shape in 2D or 2D shape in 3D, which lies on the plane
                               and is extruded upstream to obtain the inflow particles.
                               See point 2 above for more details.
+- `average_inflow_velocity=true`: If `true`, the extrapolated inflow velocity is averaged
+                                  to impose a uniform inflow profile.
+                                  When no velocity is prescribed at the inflow,
+                                  the velocity is extrapolated from the fluid domain.
+                                  Thus, turbulent flows near the inflow can lead to
+                                  anisotropic buffer-particles distribution,
+                                  resulting in a potential numerical instability.
+                                  Averaging mitigates these effects.
 
 # Examples
 ```julia
@@ -82,18 +91,20 @@ bidirectional_flow = BoundaryZone(; plane=plane_points, plane_normal, particle_s
     This is an experimental feature and may change in any future releases.
 """
 struct BoundaryZone{BT, IC, S, ZO, ZW, FD, PN}
-    initial_condition :: IC
-    spanning_set      :: S
-    zone_origin       :: ZO
-    zone_width        :: ZW
-    flow_direction    :: FD
-    plane_normal      :: PN
-    boundary_type     :: BT
+    initial_condition       :: IC
+    spanning_set            :: S
+    zone_origin             :: ZO
+    zone_width              :: ZW
+    flow_direction          :: FD
+    plane_normal            :: PN
+    boundary_type           :: BT
+    average_inflow_velocity :: Bool
 end
 
 function BoundaryZone(; plane, plane_normal, density, particle_spacing,
                       initial_condition=nothing, extrude_geometry=nothing,
-                      open_boundary_layers::Integer, boundary_type=BidirectionalFlow())
+                      open_boundary_layers::Integer, boundary_type=BidirectionalFlow(),
+                      average_inflow_velocity=true)
     if open_boundary_layers <= 0
         throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero"))
     end
@@ -120,7 +131,8 @@ function BoundaryZone(; plane, plane_normal, density, particle_spacing,
                                       boundary_type=boundary_type)
 
     return BoundaryZone(ic, spanning_set_, zone_origin, zone_width,
-                        flow_direction, plane_normal_, boundary_type)
+                        flow_direction, plane_normal_, boundary_type,
+                        average_inflow_velocity)
 end
 
 function set_up_boundary_zone(plane, plane_normal, flow_direction, density,
diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl
index b7f92787a..dd5fd1159 100644
--- a/src/schemes/boundary/open_boundary/mirroring.jl
+++ b/src/schemes/boundary/open_boundary/mirroring.jl
@@ -162,6 +162,13 @@ function extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary,
         end
     end
 
+    if !(prescribed_velocity) && boundary_zone.average_inflow_velocity
+        # When no velocity is prescribed at the inflow, the velocity is extrapolated from the fluid domain.
+        # Thus, turbulent flows near the inflow can lead to non-uniform buffer-particles distribution,
+        # resulting in a potential numerical instability. Averaging mitigates these effects.
+        average_velocity!(v_open_boundary, u_open_boundary, system, boundary_zone, semi)
+    end
+
     return system
 end
 
@@ -220,6 +227,36 @@ function mirror_position(particle_coords, boundary_zone)
     return particle_coords - 2 * dist * boundary_zone.plane_normal
 end
 
+average_velocity!(v, u, system, boundary_zone, semi) = v
+
+function average_velocity!(v, u, system, boundary_zone::BoundaryZone{InFlow}, semi)
+    (; plane_normal, zone_origin, initial_condition) = boundary_zone
+
+    # We only use the extrapolated velocity in the vicinity of the transition region.
+    # Otherwise, if the boundary zone is too large, averaging would be excessively influenced
+    # by the fluid velocity farther away from the boundary.
+    max_dist = initial_condition.particle_spacing
+
+    # This function is executed at every stage, so it is possible for buffer particles to temporarily leave the boundary zone.
+    # Thus, we use `abs()` because buffer particles may be located outside the boundary zone.
+    candidates = findall(x -> abs(dot(x - zone_origin, -plane_normal)) <= max_dist,
+                         reinterpret(reshape, SVector{ndims(system), eltype(u)},
+                                     active_coordinates(u, system)))
+
+    avg_velocity = sum(candidates) do particle
+        return current_velocity(v, system, particle) / length(candidates)
+    end
+
+    @threaded semi for particle in each_moving_particle(system)
+        # Set the velocity of the ghost node to the average velocity of the fluid domain
+        @inbounds for dim in eachindex(avg_velocity)
+            v[dim, particle] = avg_velocity[dim]
+        end
+    end
+
+    return v
+end
+
 project_velocity_on_plane_normal(vel, boundary_zone) = vel
 
 function project_velocity_on_plane_normal(vel, boundary_zone::BoundaryZone{InFlow})
diff --git a/test/schemes/boundary/open_boundary/mirroring.jl b/test/schemes/boundary/open_boundary/mirroring.jl
index 7892db837..9a81f2234 100644
--- a/test/schemes/boundary/open_boundary/mirroring.jl
+++ b/test/schemes/boundary/open_boundary/mirroring.jl
@@ -57,6 +57,7 @@
         @testset verbose=true "plane normal $i" for i in eachindex(files)
             inflow = BoundaryZone(; plane=plane_boundary[i], boundary_type=InFlow(),
                                   plane_normal=plane_boundary_normal[i],
+                                  average_inflow_velocity=false,
                                   open_boundary_layers=10, density=1000.0, particle_spacing)
 
             open_boundary = OpenBoundarySPHSystem(inflow; fluid_system,
@@ -152,6 +153,7 @@
         @testset verbose=true "plane normal $i" for i in eachindex(files)
             inflow = BoundaryZone(; plane=plane_boundary[i], boundary_type=InFlow(),
                                   plane_normal=plane_boundary_normal[i],
+                                  average_inflow_velocity=false,
                                   open_boundary_layers=10, density=1000.0, particle_spacing)
 
             open_boundary = OpenBoundarySPHSystem(inflow; fluid_system,
@@ -190,4 +192,63 @@
             @test isapprox(open_boundary.pressure, expected_pressure, atol=1e-2)
         end
     end
+
+    @testset verbose=true "Average Inflow Velocity $i-D" for i in (2, 3)
+        particle_spacing = 0.05
+        domain_length = 1.0
+        open_boundary_layers = 40
+
+        n_particles_xy = round(Int, domain_length / particle_spacing)
+
+        if i == 2
+            domain_fluid = RectangularShape(particle_spacing, (2, 1) .* n_particles_xy,
+                                            (0.0, 0.0), density=1000.0,
+                                            velocity=x -> SVector{2}(x[1], 0.0))
+
+        else
+            domain_fluid = RectangularShape(particle_spacing, (2, 1, 1) .* n_particles_xy,
+                                            (0.0, 0.0, 0.0), density=1000.0,
+                                            velocity=x -> SVector{3}(x[1], 0.0, 0.0))
+        end
+
+        smoothing_length = 1.5 * particle_spacing
+        smoothing_kernel = WendlandC2Kernel{ndims(domain_fluid)}()
+        fluid_system = EntropicallyDampedSPHSystem(domain_fluid, smoothing_kernel,
+                                                   smoothing_length, 1.0)
+        fluid_system.cache.density .= 1000.0
+
+        if i == 2
+            plane_in = ([0.0, 0.0], [0.0, domain_length])
+        else
+            plane_in = ([0.0, 0.0, 0.0], [0.0, domain_length, 0.0],
+                        [0.0, 0.0, domain_length])
+        end
+
+        inflow = BoundaryZone(; plane=plane_in, boundary_type=InFlow(),
+                              plane_normal=(i == 2 ? [1.0, 0.0] : [1.0, 0.0, 0.0]),
+                              open_boundary_layers=open_boundary_layers, density=1000.0,
+                              particle_spacing, average_inflow_velocity=true)
+        open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system,
+                                                 boundary_model=BoundaryModelTafuni(),
+                                                 buffer_size=0)
+
+        semi = Semidiscretization(fluid_system, open_boundary_in)
+        TrixiParticles.initialize_neighborhood_searches!(semi)
+
+        v_open_boundary = zero(inflow.initial_condition.velocity)
+        v_fluid = vcat(domain_fluid.velocity, domain_fluid.pressure')
+
+        TrixiParticles.set_zero!(open_boundary_in.pressure)
+
+        TrixiParticles.extrapolate_values!(open_boundary_in, v_open_boundary, v_fluid,
+                                           inflow.initial_condition.coordinates,
+                                           domain_fluid.coordinates, semi, 0.0)
+
+        # Since the velocity profile increases linearly in positive x-direction,
+        # we can use the first velocity entry as a representative value.
+        v_x_fluid_first = v_fluid[1, 1]
+
+        @test length(unique(v_open_boundary[1, :])) == 1
+        @test isapprox(-v_x_fluid_first, first(unique(v_open_boundary[1, :])))
+    end
 end