diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_viscous_shock.jl b/examples/tree_2d_dgsem/elixir_navierstokes_viscous_shock.jl
index 6a06e3cb4f6..5f34292a840 100644
--- a/examples/tree_2d_dgsem/elixir_navierstokes_viscous_shock.jl
+++ b/examples/tree_2d_dgsem/elixir_navierstokes_viscous_shock.jl
@@ -153,6 +153,7 @@ boundary_conditions_parabolic = (x_neg = boundary_condition_parabolic,
 
 semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
                                              initial_condition, solver;
+                                             solver_parabolic = ViscousFormulationBassiRebay1(),
                                              boundary_conditions = (boundary_conditions,
                                                                     boundary_conditions_parabolic))
 
diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl
index a944fdca4b6..6c3a8069a32 100644
--- a/src/solvers/dgsem_tree/dg_2d.jl
+++ b/src/solvers/dgsem_tree/dg_2d.jl
@@ -536,43 +536,43 @@ function prolong2interfaces!(cache, u, mesh::TreeMesh{2}, equations,
 
         if orientations[interface] == 1
             # interface in x-direction
-            for j in eachnode(dg), v in eachvariable(equations)
-                # Interpolate to the interfaces using a local variable for
-                # the accumulation of values (to reduce global memory operations).
-                interface_u_1 = zero(eltype(interfaces_u))
-                interface_u_2 = zero(eltype(interfaces_u))
-                for ii in eachnode(dg)
-                    # Not += to allow `@muladd` to turn these into FMAs
-                    # (see comment at the top of the file)
-                    interface_u_1 = (interface_u_1 +
-                                     u[v, ii, j, left_element] *
-                                     boundary_interpolation[ii, 2])
-                    interface_u_2 = (interface_u_2 +
-                                     u[v, ii, j, right_element] *
-                                     boundary_interpolation[ii, 1])
+            for j in eachnode(dg)
+                for v in eachvariable(equations)
+                    # Interpolate to the interfaces using a local variable for
+                    # the accumulation of values (to reduce global memory operations).
+                    interface_u_1 = zero(eltype(interfaces_u))
+                    interface_u_2 = zero(eltype(interfaces_u))
+                    for ii in eachnode(dg)
+                        # Not += to allow `@muladd` to turn these into FMAs
+                        # (see comment at the top of the file)
+                        interface_u_1 = (interface_u_1 +
+                                         u[v, ii, j, left_element] *
+                                         boundary_interpolation[ii, 2])
+                        interface_u_2 = (interface_u_2 +
+                                         u[v, ii, j, right_element] *
+                                         boundary_interpolation[ii, 1])
+                    end
+                    interfaces_u[1, v, j, interface] = interface_u_1
+                    interfaces_u[2, v, j, interface] = interface_u_2
                 end
-                interfaces_u[1, v, j, interface] = interface_u_1
-                interfaces_u[2, v, j, interface] = interface_u_2
             end
         else # if orientations[interface] == 2
             # interface in y-direction
-            for i in eachnode(dg), v in eachvariable(equations)
-                # Interpolate to the interfaces using a local variable for
-                # the accumulation of values (to reduce global memory operations).
-                interface_u_1 = zero(eltype(interfaces_u))
-                interface_u_2 = zero(eltype(interfaces_u))
-                for jj in eachnode(dg)
-                    # Not += to allow `@muladd` to turn these into FMAs
-                    # (see comment at the top of the file)
-                    interface_u_1 = (interface_u_1 +
-                                     u[v, i, jj, left_element] *
-                                     boundary_interpolation[jj, 2])
-                    interface_u_2 = (interface_u_2 +
-                                     u[v, i, jj, right_element] *
-                                     boundary_interpolation[jj, 1])
+            for i in eachnode(dg)
+                for v in eachvariable(equations)
+                    interface_u_1 = zero(eltype(interfaces_u))
+                    interface_u_2 = zero(eltype(interfaces_u))
+                    for jj in eachnode(dg)
+                        interface_u_1 = (interface_u_1 +
+                                         u[v, i, jj, left_element] *
+                                         boundary_interpolation[jj, 2])
+                        interface_u_2 = (interface_u_2 +
+                                         u[v, i, jj, right_element] *
+                                         boundary_interpolation[jj, 1])
+                    end
+                    interfaces_u[1, v, i, interface] = interface_u_1
+                    interfaces_u[2, v, i, interface] = interface_u_2
                 end
-                interfaces_u[1, v, i, interface] = interface_u_1
-                interfaces_u[2, v, i, interface] = interface_u_2
             end
         end
     end
@@ -713,52 +713,60 @@ function prolong2boundaries!(cache, u,
             # boundary in x-direction
             if neighbor_sides[boundary] == 1
                 # element in -x direction of boundary => interpolate to right boundary node (+1)
-                for l in eachnode(dg), v in eachvariable(equations)
-                    # Interpolate to the boundaries using a local variable for
-                    # the accumulation of values (to reduce global memory operations).
-                    boundary_u = zero(eltype(boundaries.u))
-                    for ii in eachnode(dg)
-                        # Not += to allow `@muladd` to turn these into FMAs
-                        # (see comment at the top of the file)
-                        boundary_u = (boundary_u +
-                                      u[v, ii, l, element] *
-                                      boundary_interpolation[ii, 2])
+                for l in eachnode(dg)
+                    for v in eachvariable(equations)
+                        # Interpolate to the boundaries using a local variable for
+                        # the accumulation of values (to reduce global memory operations).
+                        boundary_u = zero(eltype(boundaries.u))
+                        for ii in eachnode(dg)
+                            # Not += to allow `@muladd` to turn these into FMAs
+                            # (see comment at the top of the file)
+                            boundary_u = (boundary_u +
+                                          u[v, ii, l, element] *
+                                          boundary_interpolation[ii, 2])
+                        end
+                        boundaries.u[1, v, l, boundary] = boundary_u
                     end
-                    boundaries.u[1, v, l, boundary] = boundary_u
                 end
             else # element in +x direction of boundary => interpolate to left boundary node (-1)
-                for l in eachnode(dg), v in eachvariable(equations)
-                    boundary_u = zero(eltype(boundaries.u))
-                    for ii in eachnode(dg)
-                        boundary_u = (boundary_u +
-                                      u[v, ii, l, element] *
-                                      boundary_interpolation[ii, 1])
+                for l in eachnode(dg)
+                    for v in eachvariable(equations)
+                        boundary_u = zero(eltype(boundaries.u))
+                        for ii in eachnode(dg)
+                            boundary_u = (boundary_u +
+                                          u[v, ii, l, element] *
+                                          boundary_interpolation[ii, 1])
+                        end
+                        boundaries.u[2, v, l, boundary] = boundary_u
                     end
-                    boundaries.u[2, v, l, boundary] = boundary_u
                 end
             end
         else # if orientations[boundary] == 2
             # boundary in y-direction
             if neighbor_sides[boundary] == 1
                 # element in -y direction of boundary => interpolate to right boundary node (+1)
-                for l in eachnode(dg), v in eachvariable(equations)
-                    boundary_u = zero(eltype(boundaries.u))
-                    for ii in eachnode(dg)
-                        boundary_u = (boundary_u +
-                                      u[v, l, ii, element] *
-                                      boundary_interpolation[ii, 2])
+                for l in eachnode(dg)
+                    for v in eachvariable(equations)
+                        boundary_u = zero(eltype(boundaries.u))
+                        for jj in eachnode(dg)
+                            boundary_u = (boundary_u +
+                                          u[v, l, jj, element] *
+                                          boundary_interpolation[jj, 2])
+                        end
+                        boundaries.u[1, v, l, boundary] = boundary_u
                     end
-                    boundaries.u[1, v, l, boundary] = boundary_u
                 end
             else # element in +y direction of boundary => interpolate to left boundary node (-1)
-                for l in eachnode(dg), v in eachvariable(equations)
-                    boundary_u = zero(eltype(boundaries.u))
-                    for ii in eachnode(dg)
-                        boundary_u = (boundary_u +
-                                      u[v, l, ii, element] *
-                                      boundary_interpolation[ii, 1])
+                for l in eachnode(dg)
+                    for v in eachvariable(equations)
+                        boundary_u = zero(eltype(boundaries.u))
+                        for jj in eachnode(dg)
+                            boundary_u = (boundary_u +
+                                          u[v, l, jj, element] *
+                                          boundary_interpolation[jj, 1])
+                        end
+                        boundaries.u[2, v, l, boundary] = boundary_u
                     end
-                    boundaries.u[2, v, l, boundary] = boundary_u
                 end
             end
         end
@@ -1113,6 +1121,16 @@ function calc_mortar_flux!(surface_flux_values,
     return nothing
 end
 
+# For Gauss-Legendre DGSEM mortars are not yet implemented
+function calc_mortar_flux!(surface_flux_values,
+                           mesh::TreeMesh{2},
+                           have_nonconservative_terms, equations,
+                           mortar::Nothing, surface_integral,
+                           dg::DGSEM{<:GaussLegendreBasis}, cache)
+    @assert isempty(eachmortar(dg, cache))
+    return nothing
+end
+
 @inline function calc_fstar!(destination::AbstractArray{<:Any, 2}, equations,
                              surface_flux, dg::DGSEM,
                              u_interfaces, interface, orientation)
diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
index d7e452a97db..42930688aff 100644
--- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
@@ -208,6 +208,8 @@ function prolong2interfaces!(cache, flux_viscous::Tuple,
                              dg::DG)
     @unpack interfaces = cache
     @unpack orientations, neighbor_ids = interfaces
+
+    # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*!
     interfaces_u = interfaces.u
 
     flux_viscous_x, flux_viscous_y = flux_viscous
@@ -219,7 +221,6 @@ function prolong2interfaces!(cache, flux_viscous::Tuple,
         if orientations[interface] == 1
             # interface in x-direction
             for j in eachnode(dg), v in eachvariable(equations_parabolic)
-                # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*!
                 interfaces_u[1, v, j, interface] = flux_viscous_x[v, nnodes(dg), j,
                                                                   left_element]
                 interfaces_u[2, v, j, interface] = flux_viscous_x[v, 1, j,
@@ -228,7 +229,6 @@ function prolong2interfaces!(cache, flux_viscous::Tuple,
         else # if orientations[interface] == 2
             # interface in y-direction
             for i in eachnode(dg), v in eachvariable(equations_parabolic)
-                # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*!
                 interfaces_u[1, v, i, interface] = flux_viscous_y[v, i, nnodes(dg),
                                                                   left_element]
                 interfaces_u[2, v, i, interface] = flux_viscous_y[v, i, 1,
@@ -240,6 +240,77 @@ function prolong2interfaces!(cache, flux_viscous::Tuple,
     return nothing
 end
 
+# This is the version used when calculating the divergence of the viscous fluxes.
+# Specialization `flux_viscous::Tuple` needed to
+# avoid amibiguity with the hyperbolic version of `prolong2interfaces!` in dg_2d.jl
+# which is for the variables itself, i.e., `u::Array{uEltype, 4}`.
+function prolong2interfaces!(cache, flux_viscous::Tuple,
+                             mesh::TreeMesh{2},
+                             equations_parabolic::AbstractEquationsParabolic,
+                             dg::DGSEM{<:GaussLegendreBasis})
+    @unpack interfaces = cache
+    @unpack orientations, neighbor_ids = interfaces
+    @unpack boundary_interpolation = dg.basis
+
+    # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*!
+    interfaces_u = interfaces.u
+
+    flux_viscous_x, flux_viscous_y = flux_viscous
+
+    @threaded for interface in eachinterface(dg, cache)
+        left_element = neighbor_ids[1, interface]
+        right_element = neighbor_ids[2, interface]
+
+        if orientations[interface] == 1
+            # interface in x-direction
+            for j in eachnode(dg)
+                for v in eachvariable(equations_parabolic)
+                    # Interpolate to the interfaces using a local variable for
+                    # the accumulation of values (to reduce global memory operations).
+                    interface_u_1 = zero(eltype(interfaces_u))
+                    interface_u_2 = zero(eltype(interfaces_u))
+                    for ii in eachnode(dg)
+                        # Not += to allow `@muladd` to turn these into FMAs
+                        # (see comment at the top of the file)
+                        # Need `boundary_interpolation` at right (+1) node for left element
+                        interface_u_1 = (interface_u_1 +
+                                         flux_viscous_x[v, ii, j, left_element] *
+                                         boundary_interpolation[ii, 2])
+                        # Need `boundary_interpolation` at left (-1) node for right element
+                        interface_u_2 = (interface_u_2 +
+                                         flux_viscous_x[v, ii, j, right_element] *
+                                         boundary_interpolation[ii, 1])
+                    end
+                    interfaces_u[1, v, j, interface] = interface_u_1
+                    interfaces_u[2, v, j, interface] = interface_u_2
+                end
+            end
+        else # if orientations[interface] == 2
+            # interface in y-direction
+            for i in eachnode(dg)
+                for v in eachvariable(equations_parabolic)
+                    interface_u_1 = zero(eltype(interfaces_u))
+                    interface_u_2 = zero(eltype(interfaces_u))
+                    for jj in eachnode(dg)
+                        # Need `boundary_interpolation` at right (+1) node for left element
+                        interface_u_1 = (interface_u_1 +
+                                         flux_viscous_y[v, i, jj, left_element] *
+                                         boundary_interpolation[jj, 2])
+                        # Need `boundary_interpolation` at left (-1) node for right element
+                        interface_u_2 = (interface_u_2 +
+                                         flux_viscous_y[v, i, jj, right_element] *
+                                         boundary_interpolation[jj, 1])
+                    end
+                    interfaces_u[1, v, i, interface] = interface_u_1
+                    interfaces_u[2, v, i, interface] = interface_u_2
+                end
+            end
+        end
+    end
+
+    return nothing
+end
+
 # This is the version used when calculating the divergence of the viscous fluxes
 function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{2},
                               equations_parabolic, dg::DG, parabolic_scheme,
@@ -287,6 +358,8 @@ function prolong2boundaries!(cache, flux_viscous::Tuple,
                              dg::DG)
     @unpack boundaries = cache
     @unpack orientations, neighbor_sides, neighbor_ids = boundaries
+
+    # OBS! `boundaries_u` stores the "interpolated" *fluxes* and *not the solution*!
     boundaries_u = boundaries.u
     flux_viscous_x, flux_viscous_y = flux_viscous
 
@@ -298,13 +371,11 @@ function prolong2boundaries!(cache, flux_viscous::Tuple,
             if neighbor_sides[boundary] == 1
                 # element in -x direction of boundary
                 for l in eachnode(dg), v in eachvariable(equations_parabolic)
-                    # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
                     boundaries_u[1, v, l, boundary] = flux_viscous_x[v, nnodes(dg), l,
                                                                      element]
                 end
             else # Element in +x direction of boundary
                 for l in eachnode(dg), v in eachvariable(equations_parabolic)
-                    # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
                     boundaries_u[2, v, l, boundary] = flux_viscous_x[v, 1, l,
                                                                      element]
                 end
@@ -314,14 +385,12 @@ function prolong2boundaries!(cache, flux_viscous::Tuple,
             if neighbor_sides[boundary] == 1
                 # element in -y direction of boundary
                 for l in eachnode(dg), v in eachvariable(equations_parabolic)
-                    # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
                     boundaries_u[1, v, l, boundary] = flux_viscous_y[v, l, nnodes(dg),
                                                                      element]
                 end
             else
                 # element in +y direction of boundary
                 for l in eachnode(dg), v in eachvariable(equations_parabolic)
-                    # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
                     boundaries_u[2, v, l, boundary] = flux_viscous_y[v, l, 1,
                                                                      element]
                 end
@@ -332,6 +401,91 @@ function prolong2boundaries!(cache, flux_viscous::Tuple,
     return nothing
 end
 
+# This is the version used when calculating the divergence of the viscous fluxes.
+# Specialization `flux_viscous::Tuple` needed to
+# avoid amibiguity with the hyperbolic version of `prolong2boundaries!` in dg_2d.jl
+# which is for the variables itself, i.e., `u::Array{uEltype, 4}`.
+function prolong2boundaries!(cache, flux_viscous::Tuple,
+                             mesh::TreeMesh{2},
+                             equations_parabolic::AbstractEquationsParabolic,
+                             dg::DGSEM{<:GaussLegendreBasis})
+    @unpack boundaries = cache
+    @unpack orientations, neighbor_sides, neighbor_ids = boundaries
+    @unpack boundary_interpolation = dg.basis
+
+    # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
+    boundaries_u = boundaries.u
+    flux_viscous_x, flux_viscous_y = flux_viscous
+
+    @threaded for boundary in eachboundary(dg, cache)
+        element = neighbor_ids[boundary]
+
+        if orientations[boundary] == 1
+            # boundary in x-direction
+            if neighbor_sides[boundary] == 1
+                # element in -x direction of boundary => interpolate to right boundary node (+1)
+                for l in eachnode(dg)
+                    for v in eachvariable(equations_parabolic)
+                        # Interpolate to the boundaries using a local variable for
+                        # the accumulation of values (to reduce global memory operations).
+                        boundary_u = zero(eltype(boundaries_u))
+                        for ii in eachnode(dg)
+                            # Not += to allow `@muladd` to turn these into FMAs
+                            # (see comment at the top of the file)
+                            boundary_u = (boundary_u +
+                                          flux_viscous_x[v, ii, l, element] *
+                                          boundary_interpolation[ii, 2])
+                        end
+                        boundaries_u[1, v, l, boundary] = boundary_u
+                    end
+                end
+            else # element in +x direction of boundary => interpolate to left boundary node (-1)
+                for l in eachnode(dg)
+                    for v in eachvariable(equations_parabolic)
+                        boundary_u = zero(eltype(boundaries_u))
+                        for ii in eachnode(dg)
+                            boundary_u = (boundary_u +
+                                          flux_viscous_x[v, ii, l, element] *
+                                          boundary_interpolation[ii, 1])
+                        end
+                        boundaries_u[2, v, l, boundary] = boundary_u
+                    end
+                end
+            end
+        else # if orientations[boundary] == 2
+            # boundary in y-direction
+            if neighbor_sides[boundary] == 1
+                # element in -y direction of boundary => interpolate to right boundary node (+1)
+                for l in eachnode(dg)
+                    for v in eachvariable(equations_parabolic)
+                        boundary_u = zero(eltype(boundaries_u))
+                        for jj in eachnode(dg)
+                            boundary_u = (boundary_u +
+                                          flux_viscous_y[v, l, jj, element] *
+                                          boundary_interpolation[jj, 2])
+                        end
+                        boundaries_u[1, v, l, boundary] = boundary_u
+                    end
+                end
+            else # element in +y direction of boundary => interpolate to left boundary node (-1)
+                for l in eachnode(dg)
+                    for v in eachvariable(equations_parabolic)
+                        boundary_u = zero(eltype(boundaries_u))
+                        for jj in eachnode(dg)
+                            boundary_u = (boundary_u +
+                                          flux_viscous_y[v, l, jj, element] *
+                                          boundary_interpolation[jj, 1])
+                        end
+                        boundaries_u[2, v, l, boundary] = boundary_u
+                    end
+                end
+            end
+        end
+    end
+
+    return nothing
+end
+
 function calc_viscous_fluxes!(flux_viscous,
                               gradients, u_transformed,
                               mesh::Union{TreeMesh{2}, P4estMesh{2}},
@@ -704,6 +858,17 @@ function calc_mortar_flux!(surface_flux_values, mesh::TreeMesh{2},
     return nothing
 end
 
+# For Gauss-Legendre DGSEM mortars are not yet implemented
+function calc_mortar_flux!(surface_flux_values, mesh::TreeMesh{2},
+                           equations_parabolic::AbstractEquationsParabolic,
+                           mortar::Nothing,
+                           surface_integral, dg::DGSEM{<:GaussLegendreBasis},
+                           parabolic_scheme, gradient_or_divergence,
+                           cache)
+    @assert isempty(eachmortar(dg, cache))
+    return nothing
+end
+
 @inline function calc_fstar!(destination::AbstractArray{<:Any, 2},
                              mesh, equations_parabolic::AbstractEquationsParabolic,
                              surface_flux, dg::DGSEM,
@@ -924,6 +1089,60 @@ function calc_surface_integral_gradient!(gradients,
     return nothing
 end
 
+function calc_surface_integral_gradient!(gradients,
+                                         mesh::TreeMesh{2}, # for dispatch only
+                                         equations_parabolic::AbstractEquationsParabolic,
+                                         dg::DGSEM{<:GaussLegendreBasis}, cache)
+    @unpack boundary_interpolation_inverse_weights = dg.basis
+    @unpack surface_flux_values = cache.elements
+
+    gradients_x, gradients_y = gradients
+
+    # Note that all fluxes have been computed with outward-pointing normal vectors.
+    # We also use explicit assignments instead of `+=` to let `@muladd` turn these
+    # into FMAs (see comment at the top of the file).
+    @threaded for element in eachelement(dg, cache)
+        for l in eachnode(dg)
+            for v in eachvariable(equations_parabolic)
+                # Aliases for repeatedly accessed variables
+                surface_flux_minus = surface_flux_values[v, l, 1, element]
+                surface_flux_plus = surface_flux_values[v, l, 2, element]
+                for ii in eachnode(dg)
+                    # surface at -x
+                    gradients_x[v, ii, l, element] = (gradients_x[v, ii, l, element] -
+                                                      surface_flux_minus *
+                                                      boundary_interpolation_inverse_weights[ii,
+                                                                                             1])
+
+                    # surface at +x
+                    gradients_x[v, ii, l, element] = (gradients_x[v, ii, l, element] +
+                                                      surface_flux_plus *
+                                                      boundary_interpolation_inverse_weights[ii,
+                                                                                             2])
+                end
+
+                surface_flux_minus = surface_flux_values[v, l, 3, element]
+                surface_flux_plus = surface_flux_values[v, l, 4, element]
+                for jj in eachnode(dg)
+                    # surface at -y
+                    gradients_y[v, l, jj, element] = (gradients_y[v, l, jj, element] -
+                                                      surface_flux_minus *
+                                                      boundary_interpolation_inverse_weights[jj,
+                                                                                             1])
+
+                    # surface at +y
+                    gradients_y[v, l, jj, element] = (gradients_y[v, l, jj, element] +
+                                                      surface_flux_plus *
+                                                      boundary_interpolation_inverse_weights[jj,
+                                                                                             2])
+                end
+            end
+        end
+    end
+
+    return nothing
+end
+
 function reset_gradients!(gradients::NTuple{2}, dg::DG, cache)
     gradients_x, gradients_y = gradients
 
diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl
index e3e763f0829..79b9aaddcec 100644
--- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl
@@ -100,6 +100,8 @@ function prolong2interfaces!(cache, flux_viscous::Tuple,
                              dg::DG)
     @unpack interfaces = cache
     @unpack orientations, neighbor_ids = interfaces
+
+    # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*!
     interfaces_u = interfaces.u
 
     flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous
@@ -112,7 +114,7 @@ function prolong2interfaces!(cache, flux_viscous::Tuple,
             # interface in x-direction
             for k in eachnode(dg), j in eachnode(dg),
                 v in eachvariable(equations_parabolic)
-                # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*!
+
                 interfaces_u[1, v, j, k, interface] = flux_viscous_x[v,
                                                                      nnodes(dg), j, k,
                                                                      left_element]
@@ -124,7 +126,7 @@ function prolong2interfaces!(cache, flux_viscous::Tuple,
             # interface in y-direction
             for k in eachnode(dg), i in eachnode(dg),
                 v in eachvariable(equations_parabolic)
-                # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*!
+
                 interfaces_u[1, v, i, k, interface] = flux_viscous_y[v,
                                                                      i, nnodes(dg), k,
                                                                      left_element]
@@ -136,7 +138,7 @@ function prolong2interfaces!(cache, flux_viscous::Tuple,
             # interface in z-direction
             for j in eachnode(dg), i in eachnode(dg),
                 v in eachvariable(equations_parabolic)
-                # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*!
+
                 interfaces_u[1, v, i, j, interface] = flux_viscous_z[v,
                                                                      i, j, nnodes(dg),
                                                                      left_element]
@@ -198,6 +200,8 @@ function prolong2boundaries!(cache, flux_viscous::Tuple,
                              dg::DG)
     @unpack boundaries = cache
     @unpack orientations, neighbor_sides, neighbor_ids = boundaries
+
+    # OBS! `boundaries_u` stores the "interpolated" *fluxes* and *not the solution*!
     boundaries_u = boundaries.u
     flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous
 
@@ -210,7 +214,7 @@ function prolong2boundaries!(cache, flux_viscous::Tuple,
                 # element in -x direction of boundary
                 for k in eachnode(dg), j in eachnode(dg),
                     v in eachvariable(equations_parabolic)
-                    # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
+
                     boundaries_u[1, v, j, k, boundary] = flux_viscous_x[v,
                                                                         nnodes(dg),
                                                                         j,
@@ -220,7 +224,7 @@ function prolong2boundaries!(cache, flux_viscous::Tuple,
             else # Element in +x direction of boundary
                 for k in eachnode(dg), j in eachnode(dg),
                     v in eachvariable(equations_parabolic)
-                    # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
+
                     boundaries_u[2, v, j, k, boundary] = flux_viscous_x[v,
                                                                         1,
                                                                         j,
@@ -234,7 +238,7 @@ function prolong2boundaries!(cache, flux_viscous::Tuple,
                 # element in -y direction of boundary
                 for k in eachnode(dg), i in eachnode(dg),
                     v in eachvariable(equations_parabolic)
-                    # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
+
                     boundaries_u[1, v, i, k, boundary] = flux_viscous_y[v,
                                                                         i,
                                                                         nnodes(dg),
@@ -245,7 +249,7 @@ function prolong2boundaries!(cache, flux_viscous::Tuple,
                 # element in +y direction of boundary
                 for k in eachnode(dg), i in eachnode(dg),
                     v in eachvariable(equations_parabolic)
-                    # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
+
                     boundaries_u[2, v, i, k, boundary] = flux_viscous_y[v,
                                                                         i,
                                                                         1,
@@ -259,7 +263,7 @@ function prolong2boundaries!(cache, flux_viscous::Tuple,
                 # element in -z direction of boundary
                 for j in eachnode(dg), i in eachnode(dg),
                     v in eachvariable(equations_parabolic)
-                    # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
+
                     boundaries_u[1, v, i, j, boundary] = flux_viscous_z[v,
                                                                         i,
                                                                         j,
@@ -270,7 +274,7 @@ function prolong2boundaries!(cache, flux_viscous::Tuple,
                 # element in +z direction of boundary
                 for j in eachnode(dg), i in eachnode(dg),
                     v in eachvariable(equations_parabolic)
-                    # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
+
                     boundaries_u[2, v, i, j, boundary] = flux_viscous_z[v,
                                                                         i,
                                                                         j,
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index dd2051cf7f4..19953579506 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -207,6 +207,19 @@ end
     @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000)
 end
 
+@trixi_testset "TreeMesh2D: elixir_advection_diffusion.jl (Gauss-Legendre)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
+                                 "elixir_advection_diffusion.jl"),
+                        solver=DGSEM(polydeg = 5, surface_flux = flux_lax_friedrichs,
+                                     basis_type = GaussLegendreBasis),
+                        initial_refinement_level=2, tspan=(0.0, 0.4),
+                        l2=[2.8254621369070895e-6], linf=[6.914648264633172e-6])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+    @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000)
+end
+
 @trixi_testset "TreeMesh2D: elixir_advection_diffusion.jl (LDG)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
                                  "elixir_advection_diffusion.jl"),
@@ -316,6 +329,20 @@ end
     @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000)
 end
 
+@trixi_testset "TreeMesh2D: elixir_advection_diffusion_nonperiodic.jl (Gauss-Legendre)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
+                                 "elixir_advection_diffusion_nonperiodic.jl"),
+                        solver=DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs,
+                                     basis_type = GaussLegendreBasis),
+                        initial_refinement_level=2, tspan=(0.0, 0.1),
+                        l2=[0.005916696880764326],
+                        linf=[0.034212013224034776])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+    @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000)
+end
+
 @trixi_testset "TreeMesh2D: elixir_advection_diffusion_nonperiodic_amr.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
                                  "elixir_advection_diffusion_nonperiodic_amr.jl"),
@@ -672,6 +699,31 @@ end
     @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000)
 end
 
+@trixi_testset "TreeMesh2D: elixir_navierstokes_viscous_shock.jl (Gauss-Legendre, LDG)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
+                                 "elixir_navierstokes_viscous_shock.jl"),
+                        solver=DGSEM(polydeg = 3, surface_flux = flux_hlle,
+                                     basis_type = GaussLegendreBasis),
+                        solver_parabolic=ViscousFormulationLocalDG(),
+                        cfl_diffusive=0.04,
+                        l2=[
+                            6.599006355897759e-6,
+                            4.514805201434994e-6,
+                            6.54834144833621e-17,
+                            4.882545625516753e-6
+                        ],
+                        linf=[
+                            3.7580718253771295e-5,
+                            2.6691756676799905e-5,
+                            3.560074538214949e-16,
+                            2.989434893274634e-5
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+    @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000)
+end
+
 @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem",
                                  "elixir_advection_diffusion_periodic.jl"),