diff --git a/NEWS.md b/NEWS.md
index 3d64f7b78f..cba0d3f86d 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -6,13 +6,27 @@ for human readability.
 
 ## Changes when updating to v0.9 from v0.8.x
 
+#### Added
+
+- Boundary conditions are now supported on nonconservative terms ([#2062]).
+
 #### Changed
 
 - We removed the first argument `semi` corresponding to a `Semidiscretization` from the 
   `AnalysisSurfaceIntegral` constructor, as it is no longer needed (see [#1959]).
-  The `AnalysisSurfaceIntegral` now only takes the arguments `boundary_symbols` and `variable`. ([#2069])
+  The `AnalysisSurfaceIntegral` now only takes the arguments `boundary_symbols` and `variable`.
+  ([#2069])
 - In functions `rhs!`, `rhs_parabolic!`  we removed the unused argument `initial_condition`. ([#2037])
   Users should not be affected by this.
+- Nonconservative terms depend only on `normal_direction_average` instead of both 
+  `normal_direction_average` and `normal_direction_ll`, such that the function signature is now 
+  identical with conservative fluxes. This required a change of the `normal_direction` in
+  `flux_nonconservative_powell` ([#2062]).
+
+#### Deprecated
+
+#### Removed
+
 
 ## Changes in the v0.8 lifecycle
 
diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl
index d3aa4ecf58..de274248b4 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/basic_types.jl b/src/basic_types.jl
index ee479a6203..8a2430c1a8 100644
--- a/src/basic_types.jl
+++ b/src/basic_types.jl
@@ -71,14 +71,14 @@ struct BoundaryConditionDoNothing end
                                                 orientation_or_normal_direction,
                                                 direction::Integer, x, t, surface_flux,
                                                 equations)
-    return flux(u_inner, orientation_or_normal_direction, equations)
+    return surface_flux(u_inner, u_inner, orientation_or_normal_direction, equations)
 end
 
 # This version can be called by hyperbolic solvers on unstructured, curved meshes
 @inline function (::BoundaryConditionDoNothing)(u_inner,
                                                 outward_direction::AbstractVector,
                                                 x, t, surface_flux, equations)
-    return flux(u_inner, outward_direction, equations)
+    return surface_flux(u_inner, u_inner, outward_direction, equations)
 end
 
 # This version can be called by parabolic solvers
diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl
index 7a6d19facd..a43c6b0e61 100644
--- a/src/equations/ideal_glm_mhd_2d.jl
+++ b/src/equations/ideal_glm_mhd_2d.jl
@@ -196,18 +196,15 @@ 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
 Powell and the Galilean nonconservative term associated with the GLM multiplier
 of the [`IdealGlmMhdEquations2D`](@ref).
 
-On curvilinear meshes, this nonconservative flux depends on both the
-contravariant vector (normal direction) at the current node and the averaged
-one. This is different from numerical fluxes used to discretize conservative
-terms.
+On curvilinear meshes, the implementation differs from the reference since we use the averaged 
+normal direction for the metrics dealiasing.
 
 ## References
 - Marvin Bohm, Andrew R.Winters, Gregor J. Gassner, Dominik Derigs,
@@ -254,8 +251,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
@@ -265,14 +261,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})
@@ -300,8 +291,9 @@ of the [`IdealGlmMhdEquations2D`](@ref).
 
 This implementation uses a non-conservative term that can be written as the product
 of local and symmetric parts. It is equivalent to the non-conservative flux of Bohm
-et al. (`flux_nonconservative_powell`) for conforming meshes but it yields different
-results on non-conforming meshes(!).
+et al. [`flux_nonconservative_powell`](@ref) for conforming meshes but it yields different
+results on non-conforming meshes(!). On curvilinear meshes this formulation applies the
+local normal direction compared to the averaged one used in [`flux_nonconservative_powell`](@ref).
 
 The two other flux functions with the same name return either the local
 or symmetric portion of the non-conservative flux based on the type of the
diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl
index e922a2e6fd..9bb05cb3f8 100644
--- a/src/equations/ideal_glm_mhd_3d.jl
+++ b/src/equations/ideal_glm_mhd_3d.jl
@@ -224,18 +224,15 @@ 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
 Powell and the Galilean nonconservative term associated with the GLM multiplier
 of the [`IdealGlmMhdEquations3D`](@ref).
 
-On curvilinear meshes, this nonconservative flux depends on both the
-contravariant vector (normal direction) at the current node and the averaged
-one. This is different from numerical fluxes used to discretize conservative
-terms.
+On curvilinear meshes, the implementation differs from the reference since we use the averaged 
+normal direction for the metrics dealiasing.
 
 ## References
 - Marvin Bohm, Andrew R.Winters, Gregor J. Gassner, Dominik Derigs,
@@ -292,8 +289,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
@@ -303,16 +299,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 4ecaf3b6e1..30eca095aa 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 be7dd3241d..2d588c5c79 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 a2edf4ceca..88f0660701 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 17b9af0446..3c86828918 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 ece4840b74..5aabbf7ac6 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 39bece53ca..cb4e75149e 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 616a6b1513..aeae63183b 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 1f8c1ffa16..48d4fe153c 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_dgmulti_2d.jl b/test/test_dgmulti_2d.jl
index 932fb9bc95..2e41591d52 100644
--- a/test/test_dgmulti_2d.jl
+++ b/test/test_dgmulti_2d.jl
@@ -788,26 +788,26 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_reflective_wall.jl"),
                         cells_per_dimension=4,
                         l2=[
-                            0.0036019536614619687,
-                            0.001734097206958611,
-                            0.008375221008997178,
+                            0.0036019562526881602,
+                            0.0017340971255535853,
+                            0.00837522167692243,
                             0.0,
-                            0.028596796602124414,
-                            0.0018573693138866614,
-                            0.0020807798141551166,
+                            0.028596802654003512,
+                            0.0018573697892233679,
+                            0.0020807798940528956,
                             0.0,
-                            5.301188920230166e-5
+                            5.301259762428258e-5
                         ],
                         linf=[
-                            0.01692601228199253,
-                            0.009369662298436778,
-                            0.04145169295835428,
+                            0.016925983823703028,
+                            0.009369659529710701,
+                            0.04145170727840005,
                             0.0,
-                            0.11569908670112738,
-                            0.00984964453299233,
-                            0.01141708032148614,
+                            0.1156990108418654,
+                            0.009849648257876749,
+                            0.011417088537145403,
                             0.0,
-                            0.0002992631411931389
+                            0.0002992621756946904
                         ])
     # Ensure that we do not have excessive memory allocations 
     # (e.g., from type instabilities) 
diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl
index cba79b58b7..46dfa542f3 100644
--- a/test/test_p4est_2d.jl
+++ b/test/test_p4est_2d.jl
@@ -617,17 +617,18 @@ end
 
 @trixi_testset "elixir_mhd_rotor.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"),
-                        l2=[0.4552094211937344, 0.8918052934760807, 0.8324715234680768,
+                        l2=[0.45520942119628094, 0.8918052934746296, 0.8324715234698744,
                             0.0,
-                            0.9801268321975978, 0.10475722739111007,
-                            0.15551326369033164,
+                            0.9801268321975156, 0.1047572273917542,
+                            0.15551326369065485,
                             0.0,
-                            2.0602990858239632e-5],
-                        linf=[10.19421969147307, 18.254409357804683, 10.031954811332596,
+                            2.0604823850230503e-5],
+                        linf=[10.194219681102396, 18.254409373691253,
+                            10.031954811617876,
                             0.0,
-                            19.646870938371492, 1.3938679692894465, 1.8725058401937984,
+                            19.646870952133547, 1.39386796733875, 1.8725058402095736,
                             0.0,
-                            0.0016201762010257296],
+                            0.001619787048808356],
                         tspan=(0.0, 0.02))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl
index 869b7554bf..55f7ab7602 100644
--- a/test/test_p4est_3d.jl
+++ b/test/test_p4est_3d.jl
@@ -571,16 +571,16 @@ end
 
 @trixi_testset "elixir_mhd_shockcapturing_amr.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shockcapturing_amr.jl"),
-                        l2=[0.006297229188299052, 0.0064363477630573936,
-                            0.007109134822960387, 0.0065295379843073945,
-                            0.02061487028361094, 0.005561406556868266,
-                            0.007570747563219415, 0.005571060186624124,
-                            3.910359570546058e-6],
-                        linf=[0.20904050617411984, 0.18630026905465372,
-                            0.23476537952044518, 0.19430178061639747,
-                            0.6858488631108304, 0.15169972134884624,
-                            0.22431157069631724, 0.16823638724229162,
-                            0.0005352202836463904],
+                        l2=[0.006297229188267704, 0.006436347763092648,
+                            0.0071091348227321095, 0.00652953798427642,
+                            0.0206148702828057, 0.005561406556411695,
+                            0.007570747563696005, 0.005571060186513173,
+                            3.888176398720913e-6],
+                        linf=[0.20904050630623572, 0.1863002690612441,
+                            0.2347653795205547, 0.19430178062881898,
+                            0.6858488630270272, 0.15169972127018583,
+                            0.22431157058134898, 0.16823638722404644,
+                            0.0005208971463830214],
                         tspan=(0.0, 0.04),
                         coverage_override=(maxiters = 6, initial_refinement_level = 1,
                                            base_level = 1, max_level = 2))
diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl
index ba2565e0dd..f78f127f43 100644
--- a/test/test_structured_2d.jl
+++ b/test/test_structured_2d.jl
@@ -911,16 +911,16 @@ end
 
 @trixi_testset "elixir_mhd_ec.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
-                        l2=[0.04937480811868297, 0.06117033019988596,
-                            0.060998028674664716, 0.03155145889799417,
-                            0.2319175391388658, 0.02476283192966346,
-                            0.024483244374818587, 0.035439957899127385,
-                            0.0016022148194667542],
-                        linf=[0.24749024430983746, 0.2990608279625713,
-                            0.3966937932860247, 0.22265033744519683,
-                            0.9757376320946505, 0.12123736788315098,
-                            0.12837436699267113, 0.17793825293524734,
-                            0.03460761690059514],
+                        l2=[0.04937478399958968, 0.0611701500558669,
+                            0.06099805934392425, 0.031551737882277144,
+                            0.23191853685798858, 0.02476297013104899,
+                            0.024482975007695532, 0.035440179203707095,
+                            0.0016002328034991635],
+                        linf=[0.24744671083295033, 0.2990591185187605,
+                            0.3968520446251412, 0.2226544553988576,
+                            0.9752669317263143, 0.12117894533967843,
+                            0.12845218263379432, 0.17795590713819576,
+                            0.0348517136607105],
                         tspan=(0.0, 0.3))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -934,16 +934,16 @@ end
 
 @trixi_testset "elixir_mhd_alfven_wave.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
-                        l2=[0.02890769490562535, 0.0062599448721613205,
-                            0.005650300017676721, 0.007334415940022972,
-                            0.00490446035599909, 0.007202284100220619,
-                            0.007003258686714405, 0.006734267830082687,
-                            0.004253003868791559],
-                        linf=[0.17517380432288565, 0.06197353710696667,
-                            0.038494840938641646, 0.05293345499813148,
-                            0.03817506476831778, 0.042847170999492534,
-                            0.03761563456810613, 0.048184237474911844,
-                            0.04114666955364693],
+                        l2=[0.028905589451357638, 0.006259570019325034,
+                            0.005649791156739933, 0.0073272570974805004,
+                            0.004890348793116962, 0.00720944138561451,
+                            0.0069984328989438115, 0.006729800315219757,
+                            0.004318314151888631],
+                        linf=[0.17528323378978317, 0.06161030852803388,
+                            0.0388335541348234, 0.052906440559080926,
+                            0.0380036034027319, 0.04291841215471082,
+                            0.03702743958268562, 0.04815794489066357,
+                            0.0433064571343779],
                         tspan=(0.0, 1.0))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -1007,16 +1007,18 @@ end
 
 @trixi_testset "elixir_mhd_ec_shockcapturing.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec_shockcapturing.jl"),
-                        l2=[0.0364192725149364, 0.0426667193422069, 0.04261673001449095,
-                            0.025884071405646924,
-                            0.16181626564020496, 0.017346518770783536,
-                            0.017291573200291104, 0.026856206495339655,
-                            0.0007443858043598808],
-                        linf=[0.25144373906033013, 0.32881947152723745,
-                            0.3053266801502693, 0.20989755319972866,
-                            0.9927517314507455, 0.1105172121361323, 0.1257708104676617,
-                            0.1628334844841588,
-                            0.02624301627479052])
+                        l2=[0.03641928087745194, 0.04266672246194787,
+                            0.042616743034675685,
+                            0.025884076832341982,
+                            0.16181640309885276, 0.017346521291731105,
+                            0.017291600359415987, 0.026856207871456043,
+                            0.0007448774124272682],
+                        linf=[0.25144155032118376, 0.3288086335996786,
+                            0.30532573631664345, 0.20990150465080706,
+                            0.9929091025128138, 0.11053858971264774,
+                            0.12578085409726314,
+                            0.16283334251103732,
+                            0.026146463886273865])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
diff --git a/test/test_structured_3d.jl b/test/test_structured_3d.jl
index ac932b9535..f4864e45d5 100644
--- a/test/test_structured_3d.jl
+++ b/test/test_structured_3d.jl
@@ -11,7 +11,7 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "structured_3d_dgsem")
 outdir = "out"
 isdir(outdir) && rm(outdir, recursive = true)
 
-@testset "Structured mesh" begin
+@testset "StructuredMesh3D" begin
 #! format: noindent
 
 @trixi_testset "elixir_advection_basic.jl" begin
@@ -240,16 +240,17 @@ end
 
 @trixi_testset "elixir_mhd_ec.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
-                        l2=[0.009082353036644902, 0.007128360240528109,
-                            0.006970330025996491, 0.006898850266874514,
-                            0.03302008823756457, 0.003203389099143526,
-                            0.003077498677885352, 0.0030740006760477624,
-                            4.192129696970217e-5],
-                        linf=[0.2883946030582689, 0.25956437344015054,
-                            0.2614364943543665, 0.24617277938134657,
-                            1.1370443512475847, 0.1278041831463388, 0.13347391885068594,
-                            0.1457563463643099,
-                            0.0021174246048172563],
+                        l2=[0.009082353008355219, 0.007128360330314966,
+                            0.0069703300260751545, 0.006898850266164216,
+                            0.033020091335659474, 0.003203389281512797,
+                            0.0030774985678369746, 0.00307400076520122,
+                            4.192572922118587e-5],
+                        linf=[0.28839460197220435, 0.25956437090703427,
+                            0.26143649456148177, 0.24617277684934058,
+                            1.1370439348603143, 0.12780410700666367,
+                            0.13347392283166903,
+                            0.145756208548534,
+                            0.0021181795153149053],
                         tspan=(0.0, 0.25))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -263,16 +264,16 @@ end
 
 @trixi_testset "elixir_mhd_alfven_wave.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
-                        l2=[0.003015476175153681, 0.00145499403283373,
-                            0.0009125744757935803, 0.0017703080480578979,
-                            0.0013046447673965966, 0.0014564863387645508,
-                            0.0013332311430907598, 0.001647832598455728,
-                            0.0013647609788548722],
-                        linf=[0.027510637768610846, 0.02797062834945721,
-                            0.01274249949295704, 0.038940694415543736,
-                            0.02200825678588325, 0.03167600959583505,
-                            0.021420957993862344, 0.03386589835999665,
-                            0.01888303191983353],
+                        l2=[0.003015390232128414, 0.0014538563096541798,
+                            0.000912478356719486, 0.0017715065044433436,
+                            0.0013017575272262197, 0.0014545437537522726,
+                            0.0013322897333898482, 0.0016493009787844212,
+                            0.0013747547738038235],
+                        linf=[0.027577067632765795, 0.027912829563483885,
+                            0.01282206030593043, 0.03911437990598213,
+                            0.021962225923304324, 0.03169774571258743,
+                            0.021591564663781426, 0.034028148178115364,
+                            0.020084593242858988],
                         # Use same polydeg as everything else to prevent long compile times in CI
                         coverage_override=(polydeg = 3,))
     # Ensure that we do not have excessive memory allocations
@@ -287,16 +288,16 @@ end
 
 @trixi_testset "elixir_mhd_alfven_wave.jl with flux_lax_friedrichs" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
-                        l2=[0.003047854479955232, 0.0014572199588782184,
-                            0.0009093737183251411, 0.0017937548694553895,
-                            0.0013010437110755424, 0.0014545607744895874,
-                            0.001328514015121245, 0.001671342529206066,
-                            0.0013653963058149186],
-                        linf=[0.027719103797310463, 0.027570111789910784,
-                            0.012561901006903103, 0.03903568568480584,
-                            0.021311996934554767, 0.03154849824135775,
-                            0.020996033645485412, 0.03403185137382961,
-                            0.019488952445771597],
+                        l2=[0.0030477691235949685, 0.00145609137038748,
+                            0.0009092809766088607, 0.0017949926915475929,
+                            0.0012981612165627713, 0.0014525841626158234,
+                            0.0013275465154956557, 0.0016728767532610933,
+                            0.0013751925705271012],
+                        linf=[0.02778552932540901, 0.027511633996169835,
+                            0.012637649797178449, 0.03920805095546112,
+                            0.02126543791857216, 0.031563506812970266,
+                            0.02116105422516923, 0.03419432640106229,
+                            0.020324891223351533],
                         surface_flux=(flux_lax_friedrichs, flux_nonconservative_powell),
                         # Use same polydeg as everything else to prevent long compile times in CI
                         coverage_override=(polydeg = 3,))
@@ -312,16 +313,16 @@ end
 
 @trixi_testset "elixir_mhd_ec_shockcapturing.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec_shockcapturing.jl"),
-                        l2=[0.009352631220872144, 0.008058649103542618,
-                            0.008027041293333663, 0.008071417851552725,
-                            0.034909149665869485, 0.00393019428600812,
-                            0.0039219074393817, 0.003906321245184237,
-                            4.197255300781248e-5],
-                        linf=[0.30749098250807516, 0.2679008863509767,
-                            0.271243087484388, 0.26545396569129537,
-                            0.9620950892188596, 0.18163281157498123,
-                            0.15995708312378454, 0.17918221526906408,
-                            0.015138346608166353],
+                        l2=[0.009352631216098996, 0.008058649096024162,
+                            0.00802704129788766, 0.008071417834885589,
+                            0.03490914976431044, 0.003930194255268652,
+                            0.003921907459117296, 0.003906321239858786,
+                            4.1971260184918575e-5],
+                        linf=[0.307491045404509, 0.26790087991041506,
+                            0.2712430701672931, 0.2654540237991884,
+                            0.9620943261873176, 0.181632512204141,
+                            0.15995711137712265, 0.1791807940466812,
+                            0.015138421396338456],
                         tspan=(0.0, 0.25),
                         # Use same polydeg as everything else to prevent long compile times in CI
                         coverage_override=(polydeg = 3,))
diff --git a/test/test_t8code_2d.jl b/test/test_t8code_2d.jl
index 644995778d..ed9edbee3d 100644
--- a/test/test_t8code_2d.jl
+++ b/test/test_t8code_2d.jl
@@ -305,16 +305,17 @@ end
     # This test is identical to the one in `test_p4est_2d.jl` besides minor
     # deviations in the expected error norms.
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"),
-                        l2=[0.44207324634847545, 0.8804644301177857, 0.8262542320669426,
+                        l2=[0.4420732463420727, 0.8804644301158163, 0.8262542320734158,
                             0.0,
-                            0.9615023124189027, 0.10386709616755131, 0.1540308191628843,
+                            0.9615023124248694, 0.10386709616933161,
+                            0.15403081916109138,
                             0.0,
-                            2.8350276854372125e-5],
-                        linf=[10.04548675437385, 17.998696852394836, 9.575802136190026,
+                            2.835066224683485e-5],
+                        linf=[10.045486750338348, 17.998696851793447, 9.57580213608948,
                             0.0,
-                            19.431290746184473, 1.3821685018474321, 1.8186235976551453,
+                            19.431290734386764, 1.3821685025605288, 1.8186235976086789,
                             0.0,
-                            0.002309422702635547],
+                            0.0023118793481168537],
                         tspan=(0.0, 0.02))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
diff --git a/test/test_threaded.jl b/test/test_threaded.jl
index 760c5ca0d7..a3d52c1923 100644
--- a/test/test_threaded.jl
+++ b/test/test_threaded.jl
@@ -227,16 +227,16 @@ end
     @trixi_testset "elixir_mhd_ec.jl" begin
         @test_trixi_include(joinpath(examples_dir(), "structured_2d_dgsem",
                                      "elixir_mhd_ec.jl"),
-                            l2=[0.04937480811868297, 0.06117033019988596,
-                                0.060998028674664716, 0.03155145889799417,
-                                0.2319175391388658, 0.02476283192966346,
-                                0.024483244374818587, 0.035439957899127385,
-                                0.0016022148194667542],
-                            linf=[0.24749024430983746, 0.2990608279625713,
-                                0.3966937932860247, 0.22265033744519683,
-                                0.9757376320946505, 0.12123736788315098,
-                                0.12837436699267113, 0.17793825293524734,
-                                0.03460761690059514],
+                            l2=[0.04937478399958968, 0.0611701500558669,
+                                0.06099805934392425, 0.031551737882277144,
+                                0.23191853685798858, 0.02476297013104899,
+                                0.024482975007695532, 0.035440179203707095,
+                                0.0016002328034991635],
+                            linf=[0.24744671083295033, 0.2990591185187605,
+                                0.3968520446251412, 0.2226544553988576,
+                                0.9752669317263143, 0.12117894533967843,
+                                0.12845218263379432, 0.17795590713819576,
+                                0.0348517136607105],
                             tspan=(0.0, 0.3))
 
         # Ensure that we do not have excessive memory allocations
diff --git a/test/test_type.jl b/test/test_type.jl
index d507cded95..62569efe0e 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
diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl
index 43b18b48e7..c433338a23 100644
--- a/test/test_unstructured_2d.jl
+++ b/test/test_unstructured_2d.jl
@@ -254,16 +254,17 @@ end
 
 @trixi_testset "elixir_mhd_ec.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
-                        l2=[0.06418293357851637, 0.12085176618704108,
-                            0.12085099342419513, 0.07743005602933221,
-                            0.1622218916638482, 0.04044434425257972,
-                            0.04044440614962498, 0.05735896706356321,
-                            0.0020992340041681734],
-                        linf=[0.1417000509328017, 0.3210578460652491, 0.335041095545175,
-                            0.22500796423572675,
-                            0.44230628074326406, 0.16743171716317784,
-                            0.16745989278866702, 0.17700588224362557,
-                            0.02692320090677309],
+                        l2=[0.06418288595515664, 0.12085170757294698,
+                            0.12085093463857763, 0.077430018507123,
+                            0.16221988122574071, 0.040444455755985195,
+                            0.04044451621612787, 0.05735903066057611,
+                            0.002095549716217215],
+                        linf=[0.14169585310190325, 0.32104342885987625,
+                            0.33503526151419405,
+                            0.22499513309636543,
+                            0.44231595436029814, 0.16750863202541477,
+                            0.1675356630213226, 0.1770099359044508,
+                            0.026783792841168948],
                         tspan=(0.0, 0.5))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -277,16 +278,16 @@ end
 
 @trixi_testset "elixir_mhd_alfven_wave.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
-                        l2=[5.377518922553881e-5, 0.09999999206243514,
-                            0.09999999206243441, 0.1414213538550799,
-                            8.770450430886394e-6, 0.0999999926130084,
-                            0.0999999926130088, 0.14142135396487032,
-                            1.1553833987291942e-5],
-                        linf=[0.00039334982566352483, 0.14144904937275282,
-                            0.14144904937277897, 0.20003315928443416,
-                            6.826863293230012e-5, 0.14146512909995967,
-                            0.14146512909994702, 0.20006706837452526,
-                            0.00013645610312810813],
+                        l2=[5.376431895349634e-5, 0.09999999205016862,
+                            0.09999999205016788, 0.14142135386740418,
+                            8.767116801867206e-6, 0.09999999259645777,
+                            0.09999999259645763, 0.14142135397626523,
+                            1.1559626795684309e-5],
+                        linf=[0.00039380173293024345, 0.14144879547840894,
+                            0.14144879547843608, 0.2000330663752416,
+                            7.021503828519293e-5, 0.14146450834000124,
+                            0.1414645083399998, 0.20006708807562765,
+                            0.0001375806459241173],
                         tspan=(0.0, 0.5))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)