diff --git a/src/lbmpy/boundaries/boundaryconditions.py b/src/lbmpy/boundaries/boundaryconditions.py
index 6ef437bf4d098edd35fd1ab0619e4af081414c71..e49edca1c980b90f592b041100ab1a20cc4fe633 100644
--- a/src/lbmpy/boundaries/boundaryconditions.py
+++ b/src/lbmpy/boundaries/boundaryconditions.py
@@ -982,10 +982,10 @@ class iMEMBounceBack(LbBoundary):
         # shifted_vel_eqs = shifted_vel_eqs.new_without_subexpressions()
         # velocity = [eq.rhs for eq in shifted_vel_eqs.new_filtered(cqc.velocity_symbols).main_assignments]
 
-        velocity_expression = [ dx4_on_dt2 * sum( c_f_ijk[i][k] * 
-                                        (- 2 / c_s_sq * wi[i] * rho * sum( [d_i * v_i for d_i, v_i in zip(c_f_ijk[i], velocity)] for j in range(self.stencil.D)
-                                                              ) / (1+q)) for i in pdfs_in_wall_dir
-                                    ) for k in range(self.stencil.D) ]
+        velocity_expression = [  dx4_on_dt2 * sum( c_f_ijk[i][k] * 
+                                    (-(2 * rho * wi[i]) /(c_s_sq * (1+q)) *
+                                        sum( [d_i * v_i for d_i, v_i in zip(c_f_ijk[i], velocity)] ) ) for i in pdfs_in_wall_dir
+                                ) for k in range(self.stencil.D) ]
                                     
         result += [ Assignment( wall_velocity[i], sp.solve(velocity_expression[i] - F_uw[i], wall_velocity[i])[0] ) for i in self.tangential_axis]
         result.append(Assignment( wall_velocity[self.normal_axis], 0 )) 
@@ -1025,8 +1025,8 @@ class iMEMBounceBack(LbBoundary):
 
         weight_info = LbmWeightInfo(lb_method, data_type=self.data_type)
         weight_of_direction = weight_info.weight_of_direction
-        f_ijk_uw = 2 / c_s_sq * rho * sum([d_i * v_i for d_i, v_i in zip(neighbor_offset, velocity)]) * weight_of_direction(
-            dir_symbol, lb_method)
+        f_ijk_uw = 2 / (c_s_sq*(1+q)) * rho * weight_of_direction( dir_symbol, lb_method) * 
+                                sum([d_i * v_i for d_i, v_i in zip(neighbor_offset, velocity)]) 
 
     # Step 6: Return Assignment
         # Eq. 20 in Asmuth2021: Note, eq. 20 should not have the 2x on LHS and should be - not +