diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
index d3aa4ecf582..de274248b45 100644
--- a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
+++ b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
@@ -31,14 +31,6 @@ equations = IdealGlmMhdEquations2D(gamma)
 
 cells_per_dimension = (32, 64)
 
-# Extend the definition of the non-conservative Powell flux functions.
-import Trixi.flux_nonconservative_powell
-function flux_nonconservative_powell(u_ll, u_rr,
-                                     normal_direction_ll::AbstractVector,
-                                     equations::IdealGlmMhdEquations2D)
-    flux_nonconservative_powell(u_ll, u_rr, normal_direction_ll, normal_direction_ll,
-                                equations)
-end
 volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
 solver = DGSEM(polydeg = 3,
                surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell),
diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl
index ab2a4b066a1..5fa2414bd66 100644
--- a/src/equations/ideal_glm_mhd_2d.jl
+++ b/src/equations/ideal_glm_mhd_2d.jl
@@ -191,8 +191,7 @@ end
     flux_nonconservative_powell(u_ll, u_rr, orientation::Integer,
                                 equations::IdealGlmMhdEquations2D)
     flux_nonconservative_powell(u_ll, u_rr,
-                                normal_direction_ll     ::AbstractVector,
-                                normal_direction_average::AbstractVector,
+                                normal_direction::AbstractVector,
                                 equations::IdealGlmMhdEquations2D)
 
 Non-symmetric two-point flux discretizing the nonconservative (source) term of
@@ -249,8 +248,7 @@ terms.
 end
 
 @inline function flux_nonconservative_powell(u_ll, u_rr,
-                                             normal_direction_ll::AbstractVector,
-                                             normal_direction_average::AbstractVector,
+                                             normal_direction::AbstractVector,
                                              equations::IdealGlmMhdEquations2D)
     rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll
     rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr
@@ -260,14 +258,9 @@ end
     v3_ll = rho_v3_ll / rho_ll
     v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
 
-    # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector
-    # at the same node location) while `B_dot_n_rr` uses the averaged normal
-    # direction. The reason for this is that `v_dot_n_ll` depends only on the left
-    # state and multiplies some gradient while `B_dot_n_rr` is used to compute
-    # the divergence of B.
-    v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2]
-    B_dot_n_rr = B1_rr * normal_direction_average[1] +
-                 B2_rr * normal_direction_average[2]
+    v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
+    B_dot_n_rr = B1_rr * normal_direction[1] +
+                 B2_rr * normal_direction[2]
 
     # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
     # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2})
diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl
index 2ffaa575243..52baf90c6a4 100644
--- a/src/equations/ideal_glm_mhd_3d.jl
+++ b/src/equations/ideal_glm_mhd_3d.jl
@@ -219,8 +219,7 @@ end
     flux_nonconservative_powell(u_ll, u_rr, orientation::Integer,
                                 equations::IdealGlmMhdEquations3D)
     flux_nonconservative_powell(u_ll, u_rr,
-                                normal_direction_ll     ::AbstractVector,
-                                normal_direction_average::AbstractVector,
+                                normal_direction::AbstractVector,
                                 equations::IdealGlmMhdEquations3D)
 
 Non-symmetric two-point flux discretizing the nonconservative (source) term of
@@ -287,8 +286,7 @@ terms.
 end
 
 @inline function flux_nonconservative_powell(u_ll, u_rr,
-                                             normal_direction_ll::AbstractVector,
-                                             normal_direction_average::AbstractVector,
+                                             normal_direction::AbstractVector,
                                              equations::IdealGlmMhdEquations3D)
     rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll
     rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr
@@ -298,16 +296,11 @@ end
     v3_ll = rho_v3_ll / rho_ll
     v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
 
-    # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector
-    # at the same node location) while `B_dot_n_rr` uses the averaged normal
-    # direction. The reason for this is that `v_dot_n_ll` depends only on the left
-    # state and multiplies some gradient while `B_dot_n_rr` is used to compute
-    # the divergence of B.
-    v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2] +
-                 v3_ll * normal_direction_ll[3]
-    B_dot_n_rr = B1_rr * normal_direction_average[1] +
-                 B2_rr * normal_direction_average[2] +
-                 B3_rr * normal_direction_average[3]
+    v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] +
+                 v3_ll * normal_direction[3]
+    B_dot_n_rr = B1_rr * normal_direction[1] +
+                 B2_rr * normal_direction[2] +
+                 B3_rr * normal_direction[3]
 
     # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
     # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2,3}, 0, 0, 0, v_{1,2,3})
diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl
index 4ecaf3b6e14..30eca095aa4 100644
--- a/src/equations/shallow_water_2d.jl
+++ b/src/equations/shallow_water_2d.jl
@@ -267,8 +267,7 @@ end
     flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer,
                                           equations::ShallowWaterEquations2D)
     flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
-                                          normal_direction_ll     ::AbstractVector,
-                                          normal_direction_average::AbstractVector,
+                                          normal_direction::AbstractVector,
                                           equations::ShallowWaterEquations2D)
 
 Non-symmetric two-point volume flux discretizing the nonconservative (source) term
@@ -303,8 +302,7 @@ Further details are available in the papers:
 end
 
 @inline function flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
-                                                       normal_direction_ll::AbstractVector,
-                                                       normal_direction_average::AbstractVector,
+                                                       normal_direction::AbstractVector,
                                                        equations::ShallowWaterEquations2D)
     # Pull the necessary left and right state information
     h_ll = waterheight(u_ll, equations)
@@ -312,8 +310,8 @@ end
 
     # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
     return SVector(0,
-                   normal_direction_average[1] * equations.gravity * h_ll * b_jump,
-                   normal_direction_average[2] * equations.gravity * h_ll * b_jump,
+                   normal_direction[1] * equations.gravity * h_ll * b_jump,
+                   normal_direction[2] * equations.gravity * h_ll * b_jump,
                    0)
 end
 
@@ -321,8 +319,7 @@ end
     flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer,
                                         equations::ShallowWaterEquations2D)
     flux_nonconservative_fjordholm_etal(u_ll, u_rr,
-                                        normal_direction_ll     ::AbstractVector,
-                                        normal_direction_average::AbstractVector,
+                                        normal_direction::AbstractVector,
                                         equations::ShallowWaterEquations2D)
 
 Non-symmetric two-point surface flux discretizing the nonconservative (source) term of
@@ -365,8 +362,7 @@ and for curvilinear 2D case in the paper:
 end
 
 @inline function flux_nonconservative_fjordholm_etal(u_ll, u_rr,
-                                                     normal_direction_ll::AbstractVector,
-                                                     normal_direction_average::AbstractVector,
+                                                     normal_direction::AbstractVector,
                                                      equations::ShallowWaterEquations2D)
     # Pull the necessary left and right state information
     h_ll, _, _, b_ll = u_ll
@@ -376,8 +372,8 @@ end
     b_jump = b_rr - b_ll
 
     # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
-    f2 = normal_direction_average[1] * equations.gravity * h_average * b_jump
-    f3 = normal_direction_average[2] * equations.gravity * h_average * b_jump
+    f2 = normal_direction[1] * equations.gravity * h_average * b_jump
+    f3 = normal_direction[2] * equations.gravity * h_average * b_jump
 
     # First and last equations do not have a nonconservative flux
     f1 = f4 = 0
@@ -424,8 +420,7 @@ end
     flux_nonconservative_audusse_etal(u_ll, u_rr, orientation::Integer,
                                       equations::ShallowWaterEquations2D)
     flux_nonconservative_audusse_etal(u_ll, u_rr,
-                                      normal_direction_ll     ::AbstractVector,
-                                      normal_direction_average::AbstractVector,
+                                      normal_direction::AbstractVector,
                                       equations::ShallowWaterEquations2D)
 
 Non-symmetric two-point surface flux that discretizes the nonconservative (source) term.
@@ -467,8 +462,7 @@ Further details for the hydrostatic reconstruction and its motivation can be fou
 end
 
 @inline function flux_nonconservative_audusse_etal(u_ll, u_rr,
-                                                   normal_direction_ll::AbstractVector,
-                                                   normal_direction_average::AbstractVector,
+                                                   normal_direction::AbstractVector,
                                                    equations::ShallowWaterEquations2D)
     # Pull the water height and bottom topography on the left
     h_ll, _, _, b_ll = u_ll
@@ -479,8 +473,8 @@ end
     # Copy the reconstructed water height for easier to read code
     h_ll_star = u_ll_star[1]
 
-    f2 = normal_direction_average[1] * equations.gravity * (h_ll^2 - h_ll_star^2)
-    f3 = normal_direction_average[2] * equations.gravity * (h_ll^2 - h_ll_star^2)
+    f2 = normal_direction[1] * equations.gravity * (h_ll^2 - h_ll_star^2)
+    f3 = normal_direction[2] * equations.gravity * (h_ll^2 - h_ll_star^2)
 
     # First and last equations do not have a nonconservative flux
     f1 = f4 = 0
diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl
index 695260f4b9b..9abda0cc303 100644
--- a/src/solvers/dgmulti/dg.jl
+++ b/src/solvers/dgmulti/dg.jl
@@ -407,11 +407,7 @@ function calc_interface_flux!(cache, surface_integral::SurfaceIntegralWeakForm,
             # Two notes on the use of `flux_nonconservative`:
             # 1. In contrast to other mesh types, only one nonconservative part needs to be
             #    computed since we loop over the elements, not the unique interfaces.
-            # 2. In general, nonconservative fluxes can depend on both the contravariant
-            #    vectors (normal direction) at the current node and the averaged ones. However,
-            #    both are the same at watertight interfaces, so we pass `normal` twice.
-            nonconservative_part = flux_nonconservative(uM, uP, normal, normal,
-                                                        equations)
+            nonconservative_part = flux_nonconservative(uM, uP, normal, equations)
             # The factor 0.5 is necessary for the nonconservative fluxes based on the
             # interpretation of global SBP operators.
             flux_face_values[idM] = (conservative_part + 0.5 * nonconservative_part) *
@@ -556,14 +552,12 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
                                                         surface_flux, equations)
 
             # Compute pointwise nonconservative numerical flux at the boundary.
-            # In general, nonconservative fluxes can depend on both the contravariant
-            # vectors (normal direction) at the current node and the averaged ones.
-            # However, there is only one `face_normal` at boundaries, which we pass in twice.
-            # Note: This does not set any type of boundary condition for the nonconservative term
-            noncons_flux_at_face_node = nonconservative_flux(u_face_values[i, f],
-                                                             u_face_values[i, f],
-                                                             face_normal, face_normal,
-                                                             equations)
+            noncons_flux_at_face_node = boundary_condition(u_face_values[i, f],
+                                                           face_normal,
+                                                           face_coordinates,
+                                                           t,
+                                                           nonconservative_flux,
+                                                           equations)
 
             flux_face_values[i, f] = (cons_flux_at_face_node +
                                       0.5 * noncons_flux_at_face_node) * Jf[i, f]
diff --git a/src/solvers/dgmulti/flux_differencing.jl b/src/solvers/dgmulti/flux_differencing.jl
index 36aa50dff4e..b49a46d0849 100644
--- a/src/solvers/dgmulti/flux_differencing.jl
+++ b/src/solvers/dgmulti/flux_differencing.jl
@@ -76,10 +76,7 @@ end
         du_i = du[i]
         for j in col_ids
             u_j = u[j]
-            # The `normal_direction::AbstractVector` has to be passed in twice.
-            # This is because on curved meshes, nonconservative fluxes are
-            # evaluated using both the normal and its average at interfaces.
-            f_ij = volume_flux(u_i, u_j, normal_direction, normal_direction, equations)
+            f_ij = volume_flux(u_i, u_j, normal_direction, equations)
             du_i = du_i + 2 * A[i, j] * f_ij
         end
         du[i] = du_i
@@ -176,11 +173,8 @@ end
         for id in nzrange(A_base, i)
             A_ij = vals[id]
             j = rows[id]
-            # The `normal_direction::AbstractVector` has to be passed in twice.
-            # This is because on curved meshes, nonconservative fluxes are
-            # evaluated using both the normal and its average at interfaces.
             u_j = u[j]
-            f_ij = volume_flux(u_i, u_j, normal_direction, normal_direction, equations)
+            f_ij = volume_flux(u_i, u_j, normal_direction, equations)
             du_i = du_i + 2 * A_ij * f_ij
         end
         du[i] = du_i
diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl
index 17b9af04467..3c868289181 100644
--- a/src/solvers/dgsem_p4est/dg_2d.jl
+++ b/src/solvers/dgsem_p4est/dg_2d.jl
@@ -223,14 +223,8 @@ end
     flux_ = surface_flux(u_ll, u_rr, normal_direction, equations)
 
     # Compute both nonconservative fluxes
-    # In general, nonconservative fluxes can depend on both the contravariant
-    # vectors (normal direction) at the current node and the averaged ones.
-    # However, both are the same at watertight interfaces, so we pass the
-    # `normal_direction` twice.
-    noncons_primary = nonconservative_flux(u_ll, u_rr, normal_direction,
-                                           normal_direction, equations)
-    noncons_secondary = nonconservative_flux(u_rr, u_ll, normal_direction,
-                                             normal_direction, equations)
+    noncons_primary = nonconservative_flux(u_ll, u_rr, normal_direction, equations)
+    noncons_secondary = nonconservative_flux(u_rr, u_ll, normal_direction, equations)
 
     # Store the flux with nonconservative terms on the primary and secondary elements
     for v in eachvariable(equations)
@@ -369,9 +363,8 @@ end
     flux_ = boundary_condition(u_inner, normal_direction, x, t, surface_flux, equations)
 
     # Compute pointwise nonconservative numerical flux at the boundary.
-    # Note: This does not set any type of boundary condition for the nonconservative term
-    noncons_ = nonconservative_flux(u_inner, u_inner, normal_direction,
-                                    normal_direction, equations)
+    noncons_ = boundary_condition(u_inner, normal_direction, x, t, nonconservative_flux,
+                                  equations)
 
     # Copy flux to element storage in the correct orientation
     for v in eachvariable(equations)
@@ -554,8 +547,7 @@ end
     # The nonconservative flux is scaled by a factor of 0.5 based on
     # the interpretation of global SBP operators coupled discontinuously via
     # central fluxes/SATs
-    noncons = nonconservative_flux(u_ll, u_rr, normal_direction, normal_direction,
-                                   equations)
+    noncons = nonconservative_flux(u_ll, u_rr, normal_direction, equations)
 
     flux_plus_noncons = flux + 0.5f0 * noncons
 
diff --git a/src/solvers/dgsem_p4est/dg_3d.jl b/src/solvers/dgsem_p4est/dg_3d.jl
index ece4840b74b..5aabbf7ac60 100644
--- a/src/solvers/dgsem_p4est/dg_3d.jl
+++ b/src/solvers/dgsem_p4est/dg_3d.jl
@@ -289,14 +289,8 @@ end
     flux_ = surface_flux(u_ll, u_rr, normal_direction, equations)
 
     # Compute both nonconservative fluxes
-    # In general, nonconservative fluxes can depend on both the contravariant
-    # vectors (normal direction) at the current node and the averaged ones.
-    # However, both are the same at watertight interfaces, so we pass the
-    # `normal_direction` twice.
-    noncons_primary = nonconservative_flux(u_ll, u_rr, normal_direction,
-                                           normal_direction, equations)
-    noncons_secondary = nonconservative_flux(u_rr, u_ll, normal_direction,
-                                             normal_direction, equations)
+    noncons_primary = nonconservative_flux(u_ll, u_rr, normal_direction, equations)
+    noncons_secondary = nonconservative_flux(u_rr, u_ll, normal_direction, equations)
 
     # Store the flux with nonconservative terms on the primary and secondary elements
     for v in eachvariable(equations)
@@ -633,8 +627,7 @@ end
     # Compute nonconservative flux and add it to the flux scaled by a factor of 0.5 based on
     # the interpretation of global SBP operators coupled discontinuously via
     # central fluxes/SATs
-    noncons = nonconservative_flux(u_ll, u_rr, normal_direction, normal_direction,
-                                   equations)
+    noncons = nonconservative_flux(u_ll, u_rr, normal_direction, equations)
     flux_plus_noncons = flux + 0.5f0 * noncons
 
     # Copy to buffer
diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl
index 6023c9e6192..035a2aa5e09 100644
--- a/src/solvers/dgsem_structured/dg_2d.jl
+++ b/src/solvers/dgsem_structured/dg_2d.jl
@@ -195,7 +195,7 @@ end
                                                    element)
             Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii)
             # Compute the contravariant nonconservative flux.
-            fluxtilde1 = nonconservative_flux(u_node, u_node_ii, Ja1_node, Ja1_avg,
+            fluxtilde1 = nonconservative_flux(u_node, u_node_ii, Ja1_avg,
                                               equations)
             integral_contribution = integral_contribution +
                                     derivative_split[i, ii] * fluxtilde1
@@ -210,7 +210,7 @@ end
             Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj)
             # compute the contravariant nonconservative flux in the direction of the
             # averaged contravariant vector
-            fluxtilde2 = nonconservative_flux(u_node, u_node_jj, Ja2_node, Ja2_avg,
+            fluxtilde2 = nonconservative_flux(u_node, u_node_jj, Ja2_avg,
                                               equations)
             integral_contribution = integral_contribution +
                                     derivative_split[j, jj] * fluxtilde2
@@ -337,11 +337,11 @@ end
             # the interpretation of global SBP operators coupled discontinuously via
             # central fluxes/SATs
             ftilde1_L = ftilde1 +
-                        0.5f0 * nonconservative_flux(u_ll, u_rr, normal_direction,
-                                             normal_direction, equations)
+                        0.5f0 *
+                        nonconservative_flux(u_ll, u_rr, normal_direction, equations)
             ftilde1_R = ftilde1 +
-                        0.5f0 * nonconservative_flux(u_rr, u_ll, normal_direction,
-                                             normal_direction, equations)
+                        0.5f0 *
+                        nonconservative_flux(u_rr, u_ll, normal_direction, equations)
 
             set_node_vars!(fstar1_L, ftilde1_L, equations, dg, i, j)
             set_node_vars!(fstar1_R, ftilde1_R, equations, dg, i, j)
@@ -377,11 +377,11 @@ end
             # the interpretation of global SBP operators coupled discontinuously via
             # central fluxes/SATs
             ftilde2_L = ftilde2 +
-                        0.5f0 * nonconservative_flux(u_ll, u_rr, normal_direction,
-                                             normal_direction, equations)
+                        0.5f0 *
+                        nonconservative_flux(u_ll, u_rr, normal_direction, equations)
             ftilde2_R = ftilde2 +
-                        0.5f0 * nonconservative_flux(u_rr, u_ll, normal_direction,
-                                             normal_direction, equations)
+                        0.5f0 *
+                        nonconservative_flux(u_rr, u_ll, normal_direction, equations)
 
             set_node_vars!(fstar2_L, ftilde2_L, equations, dg, i, j)
             set_node_vars!(fstar2_R, ftilde2_R, equations, dg, i, j)
@@ -527,18 +527,12 @@ end
         flux = sign_jacobian * surface_flux(u_ll, u_rr, normal_direction, equations)
 
         # Compute both nonconservative fluxes
-        # In general, nonconservative fluxes can depend on both the contravariant
-        # vectors (normal direction) at the current node and the averaged ones.
-        # However, both are the same at watertight interfaces, so we pass the
-        # `normal_direction` twice.
         # Scale with sign_jacobian to ensure that the normal_direction matches that
         # from the flux above
         noncons_left = sign_jacobian *
-                       nonconservative_flux(u_ll, u_rr, normal_direction,
-                                            normal_direction, equations)
+                       nonconservative_flux(u_ll, u_rr, normal_direction, equations)
         noncons_right = sign_jacobian *
-                        nonconservative_flux(u_rr, u_ll, normal_direction,
-                                             normal_direction, equations)
+                        nonconservative_flux(u_rr, u_ll, normal_direction, equations)
 
         for v in eachvariable(equations)
             # Note the factor 0.5 necessary for the nonconservative fluxes based on
diff --git a/src/solvers/dgsem_structured/dg_3d.jl b/src/solvers/dgsem_structured/dg_3d.jl
index e8430eaa491..fbf6747683b 100644
--- a/src/solvers/dgsem_structured/dg_3d.jl
+++ b/src/solvers/dgsem_structured/dg_3d.jl
@@ -228,7 +228,7 @@ end
             Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii)
             # compute the contravariant nonconservative flux in the direction of the
             # averaged contravariant vector
-            fluxtilde1 = nonconservative_flux(u_node, u_node_ii, Ja1_node, Ja1_avg,
+            fluxtilde1 = nonconservative_flux(u_node, u_node_ii, Ja1_avg,
                                               equations)
             integral_contribution = integral_contribution +
                                     derivative_split[i, ii] * fluxtilde1
@@ -243,7 +243,7 @@ end
             Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj)
             # compute the contravariant nonconservative flux in the direction of the
             # averaged contravariant vector
-            fluxtilde2 = nonconservative_flux(u_node, u_node_jj, Ja2_node, Ja2_avg,
+            fluxtilde2 = nonconservative_flux(u_node, u_node_jj, Ja2_avg,
                                               equations)
             integral_contribution = integral_contribution +
                                     derivative_split[j, jj] * fluxtilde2
@@ -258,7 +258,7 @@ end
             Ja3_avg = 0.5f0 * (Ja3_node + Ja3_node_kk)
             # compute the contravariant nonconservative flux in the direction of the
             # averaged contravariant vector
-            fluxtilde3 = nonconservative_flux(u_node, u_node_kk, Ja3_node, Ja3_avg,
+            fluxtilde3 = nonconservative_flux(u_node, u_node_kk, Ja3_avg,
                                               equations)
             integral_contribution = integral_contribution +
                                     derivative_split[k, kk] * fluxtilde3
@@ -411,11 +411,11 @@ end
             # the interpretation of global SBP operators coupled discontinuously via
             # central fluxes/SATs
             ftilde_L = ftilde +
-                       0.5f0 * nonconservative_flux(u_ll, u_rr, normal_direction,
-                                            normal_direction, equations)
+                       0.5f0 *
+                       nonconservative_flux(u_ll, u_rr, normal_direction, equations)
             ftilde_R = ftilde +
-                       0.5f0 * nonconservative_flux(u_rr, u_ll, normal_direction,
-                                            normal_direction, equations)
+                       0.5f0 *
+                       nonconservative_flux(u_rr, u_ll, normal_direction, equations)
 
             set_node_vars!(fstar1_L, ftilde_L, equations, dg, i, j, k)
             set_node_vars!(fstar1_R, ftilde_R, equations, dg, i, j, k)
@@ -449,11 +449,11 @@ end
             # the interpretation of global SBP operators coupled discontinuously via
             # central fluxes/SATs
             ftilde_L = ftilde +
-                       0.5f0 * nonconservative_flux(u_ll, u_rr, normal_direction,
-                                            normal_direction, equations)
+                       0.5f0 *
+                       nonconservative_flux(u_ll, u_rr, normal_direction, equations)
             ftilde_R = ftilde +
-                       0.5f0 * nonconservative_flux(u_rr, u_ll, normal_direction,
-                                            normal_direction, equations)
+                       0.5f0 *
+                       nonconservative_flux(u_rr, u_ll, normal_direction, equations)
 
             set_node_vars!(fstar2_L, ftilde_L, equations, dg, i, j, k)
             set_node_vars!(fstar2_R, ftilde_R, equations, dg, i, j, k)
@@ -487,11 +487,11 @@ end
             # the interpretation of global SBP operators coupled discontinuously via
             # central fluxes/SATs
             ftilde_L = ftilde +
-                       0.5f0 * nonconservative_flux(u_ll, u_rr, normal_direction,
-                                            normal_direction, equations)
+                       0.5f0 *
+                       nonconservative_flux(u_ll, u_rr, normal_direction, equations)
             ftilde_R = ftilde +
-                       0.5f0 * nonconservative_flux(u_rr, u_ll, normal_direction,
-                                            normal_direction, equations)
+                       0.5f0 *
+                       nonconservative_flux(u_rr, u_ll, normal_direction, equations)
 
             set_node_vars!(fstar3_L, ftilde_L, equations, dg, i, j, k)
             set_node_vars!(fstar3_R, ftilde_R, equations, dg, i, j, k)
@@ -663,18 +663,12 @@ end
         flux = sign_jacobian * surface_flux(u_ll, u_rr, normal_direction, equations)
 
         # Compute both nonconservative fluxes
-        # In general, nonconservative fluxes can depend on both the contravariant
-        # vectors (normal direction) at the current node and the averaged ones.
-        # However, both are the same at watertight interfaces, so we pass the
-        # `normal_direction` twice.
         # Scale with sign_jacobian to ensure that the normal_direction matches that
         # from the flux above
         noncons_left = sign_jacobian *
-                       nonconservative_flux(u_ll, u_rr, normal_direction,
-                                            normal_direction, equations)
+                       nonconservative_flux(u_ll, u_rr, normal_direction, equations)
         noncons_right = sign_jacobian *
-                        nonconservative_flux(u_rr, u_ll, normal_direction,
-                                             normal_direction, equations)
+                        nonconservative_flux(u_rr, u_ll, normal_direction, equations)
 
         for v in eachvariable(equations)
             # Note the factor 0.5 necessary for the nonconservative fluxes based on
diff --git a/src/solvers/dgsem_unstructured/dg_2d.jl b/src/solvers/dgsem_unstructured/dg_2d.jl
index a531ded91cc..40832e98fd2 100644
--- a/src/solvers/dgsem_unstructured/dg_2d.jl
+++ b/src/solvers/dgsem_unstructured/dg_2d.jl
@@ -247,14 +247,10 @@ function calc_interface_flux!(surface_flux_values,
             flux = surface_flux(u_ll, u_rr, outward_direction, equations)
 
             # Compute both nonconservative fluxes
-            # In general, nonconservative fluxes can depend on both the contravariant
-            # vectors (normal direction) at the current node and the averaged ones.
-            # However, both are the same at watertight interfaces, so we pass the
-            # `outward_direction` twice.
             noncons_primary = nonconservative_flux(u_ll, u_rr, outward_direction,
-                                                   outward_direction, equations)
+                                                   equations)
             noncons_secondary = nonconservative_flux(u_rr, u_ll, outward_direction,
-                                                     outward_direction, equations)
+                                                     equations)
 
             # Copy flux to primary and secondary element storage
             # Note the sign change for the components in the secondary element!
@@ -449,13 +445,8 @@ end
     flux = boundary_condition(u_inner, outward_direction, x, t, surface_flux, equations)
 
     # Compute pointwise nonconservative numerical flux at the boundary.
-    # In general, nonconservative fluxes can depend on both the contravariant
-    # vectors (normal direction) at the current node and the averaged ones.
-    # However, both are the same at watertight interfaces, so we pass the
-    # `outward_direction` twice.
-    # Note: This does not set any type of boundary condition for the nonconservative term
-    noncons_flux = nonconservative_flux(u_inner, u_inner, outward_direction,
-                                        outward_direction, equations)
+    noncons_flux = boundary_condition(u_inner, outward_direction, x, t,
+                                      nonconservative_flux, equations)
 
     for v in eachvariable(equations)
         # Note the factor 0.5 necessary for the nonconservative fluxes based on
diff --git a/test/test_type.jl b/test/test_type.jl
index d507cded958..62569efe0ec 100644
--- a/test/test_type.jl
+++ b/test/test_type.jl
@@ -1153,8 +1153,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                              zero(RealT))
             orientations = [1, 2]
             directions = [1, 2, 3, 4]
-            normal_direction = normal_direction_ll = normal_direction_average = SVector(one(RealT),
-                                                                                        zero(RealT))
+            normal_direction = SVector(one(RealT), zero(RealT))
             nonconservative_type_local = Trixi.NonConservativeLocal()
             nonconservative_type_symmetric = Trixi.NonConservativeSymmetric()
             nonconservative_terms = [1, 2]
@@ -1167,8 +1166,7 @@ isdir(outdir) && rm(outdir, recursive = true)
 
             @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr,
-                                                               normal_direction_ll,
-                                                               normal_direction_average,
+                                                               normal_direction,
                                                                equations)) == RealT
             @test eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, normal_direction,
                                                            equations)) == RealT
@@ -1273,9 +1271,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                              zero(RealT))
             orientations = [1, 2, 3]
             directions = [1, 2, 3, 4, 5, 6]
-            normal_direction = normal_direction_ll = normal_direction_average = SVector(one(RealT),
-                                                                                        zero(RealT),
-                                                                                        zero(RealT))
+            normal_direction = SVector(one(RealT), zero(RealT), zero(RealT))
 
             @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT
             @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) ==
@@ -1285,8 +1281,7 @@ isdir(outdir) && rm(outdir, recursive = true)
 
             @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr,
-                                                               normal_direction_ll,
-                                                               normal_direction_average,
+                                                               normal_direction,
                                                                equations)) == RealT
             @test eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, normal_direction,
                                                            equations)) == RealT
@@ -2221,8 +2216,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                        one(RealT))
             orientations = [1, 2]
             directions = [1, 2, 3, 4]
-            normal_direction = normal_direction_ll = normal_direction_average = SVector(one(RealT),
-                                                                                        zero(RealT))
+            normal_direction = SVector(one(RealT), zero(RealT))
 
             surface_flux_function = flux_lax_friedrichs
             dissipation = DissipationLocalLaxFriedrichs()
@@ -2242,17 +2236,14 @@ isdir(outdir) && rm(outdir, recursive = true)
 
             @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
             @test eltype(@inferred flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
-                                                                         normal_direction_ll,
-                                                                         normal_direction_average,
+                                                                         normal_direction,
                                                                          equations)) ==
                   RealT
             @test eltype(@inferred flux_nonconservative_fjordholm_etal(u_ll, u_rr,
-                                                                       normal_direction_ll,
-                                                                       normal_direction_average,
+                                                                       normal_direction,
                                                                        equations)) == RealT
             @test eltype(@inferred flux_nonconservative_audusse_etal(u_ll, u_rr,
-                                                                     normal_direction_ll,
-                                                                     normal_direction_average,
+                                                                     normal_direction,
                                                                      equations)) == RealT
             @test eltype(@inferred flux_fjordholm_etal(u_ll, u_rr, normal_direction,
                                                        equations)) == RealT