diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_cubed_sphere.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_cubed_sphere.jl
new file mode 100644
index 00000000000..2f9c5e09f22
--- /dev/null
+++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_cubed_sphere.jl
@@ -0,0 +1,79 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (0.0, 0.0, 0.0)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+
+########
+# For diffusion
+diffusivity() = 5.0e-2
+equations_parabolic = LaplaceDiffusion3D(diffusivity(), equations)
+########
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+polydeg = 4
+solver = DGSEM(polydeg=polydeg, surface_flux=flux_lax_friedrichs)
+
+initial_condition = initial_condition_convergence_test # A sinusoidal wave gets diffused
+#initial_condition = initial_condition_constant # A constant solution stays constant
+
+mesh = Trixi.P4estMeshCubedSphere2D(5, 1.0, polydeg=polydeg, initial_refinement_level=0)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver)
+
+# compute area of the sphere
+area = zero(Float64)
+for element in 1:Trixi.ncells(mesh)
+    for j in 1:polydeg+1, i in 1:polydeg+1
+        global area += solver.basis.weights[i] * solver.basis.weights[j] / semi.cache.elements.inverse_jacobian[i, j, element]
+    end
+end
+println("Area of sphere: ", area)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval = 10,
+                                     save_analysis = true,
+                                     extra_analysis_integrals = (Trixi.energy_total,))
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=1,
+                                     solution_variables=cons2prim)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+#stepsize_callback = StepsizeCallback(cfl=0.001)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution) #, stepsize_callback
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+#= sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0e-3, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks); =#
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+alg = RDPK3SpFSAL49()
+time_int_tol = 1.0e-11
+sol = solve(ode, alg; abstol = time_int_tol, reltol = time_int_tol,
+            ode_default_options()..., callback = callbacks)
+
+# Print the timer summary
+summary_callback()
\ No newline at end of file
diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl
index 0f44941390a..8c6c033309a 100644
--- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl
+++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl
@@ -68,7 +68,7 @@ struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic
                                                                                                            Cache,
                                                                                                            CacheParabolic
                                                                                                            }
-        @assert ndims(mesh) == ndims(equations)
+        #@assert ndims(mesh) == ndims(equations)
 
         # Todo: assert nvariables(equations)==nvariables(equations_parabolic)
 
diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index a7f3345168f..6784ab15eb0 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -10,9 +10,15 @@ function create_cache_parabolic(mesh::P4estMesh{2}, equations_hyperbolic::Abstra
     interfaces = init_interfaces(mesh, equations_hyperbolic, dg.basis, elements)
     boundaries = init_boundaries(mesh, equations_hyperbolic, dg.basis, elements)
 
-    viscous_container = init_viscous_container_2d(nvariables(equations_hyperbolic),
-                                                  nnodes(dg.basis), nelements(elements),
-                                                  uEltype)
+    if ndims(equations_hyperbolic) == 2
+        viscous_container = init_viscous_container_2d(nvariables(equations_hyperbolic),
+                                                    nnodes(dg.basis), nelements(elements),
+                                                    uEltype)
+    elseif ndims(equations_hyperbolic) == 3
+        viscous_container = init_viscous_container_2d_eq_3D(nvariables(equations_hyperbolic),
+                                                            nnodes(dg.basis), nelements(elements),
+                                                            uEltype)
+    end
 
     cache = (; elements, interfaces, boundaries, viscous_container)
 
@@ -122,8 +128,9 @@ function rhs_parabolic!(du, u, t, mesh::P4estMesh{2},
     return nothing
 end
 
+# Calculate gradient for 2D p4est meshes with 2D equations
 function calc_gradient!(gradients, u_transformed, t,
-                        mesh::P4estMesh{2}, equations_parabolic,
+                        mesh::P4estMesh{2}, equations_parabolic::AbstractEquationsParabolic{2},
                         boundary_conditions_parabolic, dg::DG,
                         cache, cache_parabolic)
     gradients_x, gradients_y = gradients
@@ -328,6 +335,256 @@ function calc_gradient!(gradients, u_transformed, t,
     return nothing
 end
 
+# Calculate gradient for 2D p4est meshes with 3D equations
+function calc_gradient!(gradients, u_transformed, t,
+                        mesh::P4estMesh{2}, equations_parabolic::AbstractEquationsParabolic{3},
+                        boundary_conditions_parabolic, dg::DG,
+                        cache, cache_parabolic)
+    gradients_x, gradients_y, gradients_z = gradients
+
+    # Reset du
+    @trixi_timeit timer() "reset gradients" begin
+        reset_du!(gradients_x, dg, cache)
+        reset_du!(gradients_y, dg, cache)
+        reset_du!(gradients_z, dg, cache)
+    end
+
+    # Calculate volume integral
+    @trixi_timeit timer() "volume integral" begin
+        (; derivative_dhat) = dg.basis
+        (; contravariant_vectors) = cache.elements
+
+        @threaded for element in eachelement(dg, cache)
+
+            # Calculate gradients with respect to reference coordinates in one element
+            for j in eachnode(dg), i in eachnode(dg)
+                u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, j,
+                                       element)
+
+                for ii in eachnode(dg)
+                    multiply_add_to_node_vars!(gradients_x, derivative_dhat[ii, i], u_node,
+                                               equations_parabolic, dg, ii, j, element)
+                end
+
+                for jj in eachnode(dg)
+                    multiply_add_to_node_vars!(gradients_y, derivative_dhat[jj, j], u_node,
+                                               equations_parabolic, dg, i, jj, element)
+                end
+            end
+
+            # now that the reference coordinate gradients are computed, transform them node-by-node to physical gradients
+            # using the contravariant vectors
+            for j in eachnode(dg), i in eachnode(dg)
+                Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors, i, j,
+                                                            element)
+                Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors, i, j,
+                                                            element)
+                Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, i, j,
+                                                            element)
+                
+                gradients_reference_1 = get_node_vars(gradients_x, equations_parabolic, dg,
+                                                      i, j, element)
+                gradients_reference_2 = get_node_vars(gradients_y, equations_parabolic, dg,
+                                                      i, j, element)
+
+                # note that the contravariant vectors are transposed compared with computations of flux
+                # divergences in `calc_volume_integral!`. See
+                # https://github.com/trixi-framework/Trixi.jl/pull/1490#discussion_r1213345190
+                # for a more detailed discussion.
+                gradient_x_node = Ja11 * gradients_reference_1 +
+                                  Ja21 * gradients_reference_2
+                gradient_y_node = Ja12 * gradients_reference_1 +
+                                  Ja22 * gradients_reference_2
+                gradient_z_node = Ja13 * gradients_reference_1 +
+                                  Ja23 * gradients_reference_2
+
+                set_node_vars!(gradients_x, gradient_x_node, equations_parabolic, dg, i, j,
+                               element)
+                set_node_vars!(gradients_y, gradient_y_node, equations_parabolic, dg, i, j,
+                               element)
+                set_node_vars!(gradients_z, gradient_z_node, equations_parabolic, dg, i, j,
+                               element)
+            end
+        end
+    end
+
+    # Prolong solution to interfaces. 
+    # This reuses `prolong2interfaces` for the purely hyperbolic case.
+    @trixi_timeit timer() "prolong2interfaces" begin
+        prolong2interfaces!(cache_parabolic, u_transformed, mesh,
+                            equations_parabolic, dg.surface_integral, dg)
+    end
+
+    # Calculate interface fluxes for the gradient. 
+    # This reuses `calc_interface_flux!` for the purely hyperbolic case.
+    @trixi_timeit timer() "interface flux" begin
+        calc_interface_flux!(cache_parabolic.elements.surface_flux_values,
+                             mesh, False(), # False() = no nonconservative terms
+                             equations_parabolic, dg.surface_integral, dg, cache_parabolic)
+    end
+
+    # Prolong solution to boundaries
+    @trixi_timeit timer() "prolong2boundaries" begin
+        prolong2boundaries!(cache_parabolic, u_transformed, mesh,
+                            equations_parabolic, dg.surface_integral, dg)
+    end
+
+    # Calculate boundary fluxes
+    @trixi_timeit timer() "boundary flux" begin
+        calc_boundary_flux_gradients!(cache_parabolic, t, boundary_conditions_parabolic,
+                                      mesh, equations_parabolic, dg.surface_integral, dg)
+    end
+
+    # Prolong solution to mortars. This resues the hyperbolic version of `prolong2mortars`
+    @trixi_timeit timer() "prolong2mortars" begin
+        prolong2mortars!(cache, u_transformed, mesh, equations_parabolic,
+                         dg.mortar, dg.surface_integral, dg)
+    end
+
+    # Calculate mortar fluxes. This reuses the hyperbolic version of `calc_mortar_flux`,
+    # along with a specialization on `calc_mortar_flux!(fstar, ...)` and `mortar_fluxes_to_elements!` for 
+    # AbstractEquationsParabolic. 
+    @trixi_timeit timer() "mortar flux" begin
+        calc_mortar_flux!(cache_parabolic.elements.surface_flux_values,
+                          mesh, False(), # False() = no nonconservative terms
+                          equations_parabolic,
+                          dg.mortar, dg.surface_integral, dg, cache)
+    end
+
+    # Calculate surface integrals
+    @trixi_timeit timer() "surface integral" begin
+        (; boundary_interpolation) = dg.basis
+        (; surface_flux_values) = cache_parabolic.elements
+        (; contravariant_vectors) = cache.elements
+
+        # Access the factors only once before beginning the loop to increase performance.
+        # We also use explicit assignments instead of `+=` to let `@muladd` turn these
+        # into FMAs (see comment at the top of the file).
+        factor_1 = boundary_interpolation[1, 1]
+        factor_2 = boundary_interpolation[nnodes(dg), 2]
+        @threaded for element in eachelement(dg, cache)
+            for l in eachnode(dg)
+                for v in eachvariable(equations_parabolic)
+
+                    # Compute x-component of gradients
+
+                    # surface at -x
+                    normal_direction_x, _ = get_normal_direction(1, contravariant_vectors,
+                                                                 1, l, element)
+                    gradients_x[v, 1, l, element] = (gradients_x[v, 1, l, element] +
+                                                     surface_flux_values[v, l, 1, element] *
+                                                     factor_1 * normal_direction_x)
+
+                    # surface at +x
+                    normal_direction_x, _ = get_normal_direction(2, contravariant_vectors,
+                                                                 nnodes(dg), l, element)
+                    gradients_x[v, nnodes(dg), l, element] = (gradients_x[v, nnodes(dg), l,
+                                                                          element] +
+                                                              surface_flux_values[v, l, 2,
+                                                                                  element] *
+                                                              factor_2 * normal_direction_x)
+
+                    # surface at -y
+                    normal_direction_x, _ = get_normal_direction(3, contravariant_vectors,
+                                                                 l, 1, element)
+                    gradients_x[v, l, 1, element] = (gradients_x[v, l, 1, element] +
+                                                     surface_flux_values[v, l, 3, element] *
+                                                     factor_1 * normal_direction_x)
+
+                    # surface at +y
+                    normal_direction_x, _ = get_normal_direction(4, contravariant_vectors,
+                                                                 l, nnodes(dg), element)
+                    gradients_x[v, l, nnodes(dg), element] = (gradients_x[v, l, nnodes(dg),
+                                                                          element] +
+                                                              surface_flux_values[v, l, 4,
+                                                                                  element] *
+                                                              factor_2 * normal_direction_x)
+
+                    # Compute y-component of gradients
+
+                    # surface at -x
+                    _, normal_direction_y = get_normal_direction(1, contravariant_vectors,
+                                                                 1, l, element)
+                    gradients_y[v, 1, l, element] = (gradients_y[v, 1, l, element] +
+                                                     surface_flux_values[v, l, 1, element] *
+                                                     factor_1 * normal_direction_y)
+
+                    # surface at +x
+                    _, normal_direction_y = get_normal_direction(2, contravariant_vectors,
+                                                                 nnodes(dg), l, element)
+                    gradients_y[v, nnodes(dg), l, element] = (gradients_y[v, nnodes(dg), l,
+                                                                          element] +
+                                                              surface_flux_values[v, l, 2,
+                                                                                  element] *
+                                                              factor_2 * normal_direction_y)
+
+                    # surface at -y
+                    _, normal_direction_y = get_normal_direction(3, contravariant_vectors,
+                                                                 l, 1, element)
+                    gradients_y[v, l, 1, element] = (gradients_y[v, l, 1, element] +
+                                                     surface_flux_values[v, l, 3, element] *
+                                                     factor_1 * normal_direction_y)
+
+                    # surface at +y
+                    _, normal_direction_y = get_normal_direction(4, contravariant_vectors,
+                                                                 l, nnodes(dg), element)
+                    gradients_y[v, l, nnodes(dg), element] = (gradients_y[v, l, nnodes(dg),
+                                                                          element] +
+                                                              surface_flux_values[v, l, 4,
+                                                                                  element] *
+                                                              factor_2 * normal_direction_y)
+
+                    # Compute z-component of gradients
+
+                    # surface at -x
+                    _, _, normal_direction_z = get_normal_direction(1, contravariant_vectors,
+                                                                    1, l, element)
+                    gradients_z[v, 1, l, element] = (gradients_z[v, 1, l, element] +
+                                                     surface_flux_values[v, l, 1, element] *
+                                                     factor_1 * normal_direction_z)
+
+                    # surface at +x
+                    _, _, normal_direction_z = get_normal_direction(2, contravariant_vectors,
+                                                                 nnodes(dg), l, element)
+                    gradients_z[v, nnodes(dg), l, element] = (gradients_z[v, nnodes(dg), l,
+                                                                          element] +
+                                                              surface_flux_values[v, l, 2,
+                                                                                  element] *
+                                                              factor_2 * normal_direction_z)
+
+                    # surface at -y
+                    _, _, normal_direction_z = get_normal_direction(3, contravariant_vectors,
+                                                                 l, 1, element)
+                    gradients_z[v, l, 1, element] = (gradients_z[v, l, 1, element] +
+                                                     surface_flux_values[v, l, 3, element] *
+                                                     factor_1 * normal_direction_z)
+
+                    # surface at +y
+                    _, _, normal_direction_z = get_normal_direction(4, contravariant_vectors,
+                                                                 l, nnodes(dg), element)
+                    gradients_z[v, l, nnodes(dg), element] = (gradients_z[v, l, nnodes(dg),
+                                                                          element] +
+                                                              surface_flux_values[v, l, 4,
+                                                                                  element] *
+                                                              factor_2 * normal_direction_z)
+                end
+            end
+        end
+    end
+
+    # Apply Jacobian from mapping to reference element
+    @trixi_timeit timer() "Jacobian" begin
+        apply_jacobian_parabolic!(gradients_x, mesh, equations_parabolic, dg,
+                                  cache_parabolic)
+        apply_jacobian_parabolic!(gradients_y, mesh, equations_parabolic, dg,
+                                  cache_parabolic)
+        apply_jacobian_parabolic!(gradients_z, mesh, equations_parabolic, dg,
+                                  cache_parabolic)
+    end
+
+    return nothing
+end
+
 # This version is called during `calc_gradients!` and must be specialized because the 
 # flux for the gradient is {u}, which doesn't depend on the outward normal. Thus, 
 # you don't need to scale by 2 (e.g., the scaling factor in the normals (and in the 
@@ -413,9 +670,10 @@ end
 end
 
 # This is the version used when calculating the divergence of the viscous fluxes
+# of a 2D equation on a 2D p4est mesh
 function calc_volume_integral!(du, flux_viscous,
                                mesh::P4estMesh{2},
-                               equations_parabolic::AbstractEquationsParabolic,
+                               equations_parabolic::AbstractEquationsParabolic{2},
                                dg::DGSEM, cache)
     (; derivative_dhat) = dg.basis
     (; contravariant_vectors) = cache.elements
@@ -451,10 +709,66 @@ function calc_volume_integral!(du, flux_viscous,
 end
 
 # This is the version used when calculating the divergence of the viscous fluxes
+# of a 3D equation on a 2D p4est mesh
+function calc_volume_integral!(du, flux_viscous,
+                               mesh::P4estMesh{2},
+                               equations_parabolic::AbstractEquationsParabolic{3},
+                               dg::DGSEM, cache)
+    (; derivative_dhat) = dg.basis
+    (; contravariant_vectors) = cache.elements
+    flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous
+
+    @threaded for element in eachelement(dg, cache)
+        # Calculate volume terms in one element
+        for j in eachnode(dg), i in eachnode(dg)
+            flux1 = get_node_vars(flux_viscous_x, equations_parabolic, dg, i, j, element)
+            flux2 = get_node_vars(flux_viscous_y, equations_parabolic, dg, i, j, element)
+            flux3 = get_node_vars(flux_viscous_z, equations_parabolic, dg, i, j, element)
+
+            # Compute the contravariant flux by taking the scalar product of the
+            # first contravariant vector Ja^1 and the flux vector
+            #Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
+            Ja11 = contravariant_vectors[1, 1, i, j, element]
+            Ja12 = contravariant_vectors[2, 1, i, j, element]
+            Ja13 = contravariant_vectors[3, 1, i, j, element]
+            contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 + Ja13 * flux3
+            for ii in eachnode(dg)
+                multiply_add_to_node_vars!(du, derivative_dhat[ii, i], contravariant_flux1,
+                                           equations_parabolic, dg, ii, j, element)
+            end
+
+            # Compute the contravariant flux by taking the scalar product of the
+            # second contravariant vector Ja^2 and the flux vector
+            #Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
+            Ja21 = contravariant_vectors[1, 2, i, j, element]
+            Ja22 = contravariant_vectors[2, 2, i, j, element]
+            Ja23 = contravariant_vectors[3, 2, i, j, element]
+            contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 + Ja23 * flux3
+            for jj in eachnode(dg)
+                multiply_add_to_node_vars!(du, derivative_dhat[jj, j], contravariant_flux2,
+                                           equations_parabolic, dg, i, jj, element)
+            end
+
+            #Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, i, j, element)
+            Ja31 = contravariant_vectors[1, 3, i, j, element]
+            Ja32 = contravariant_vectors[2, 3, i, j, element]
+            Ja33 = contravariant_vectors[3, 3, i, j, element]
+            for v in eachvariable(equations_parabolic)
+                du[v, i, j, element] += (Ja31 * flux1[v] + Ja32 * flux2[v] +
+                                         Ja33 * flux3[v])
+            end
+        end
+    end
+
+    return nothing
+end
+
+# This is the version used when calculating the divergence of the viscous fluxes
+# of a 2D equation on a 2D p4est mesh
 # We pass the `surface_integral` argument solely for dispatch
 function prolong2interfaces!(cache_parabolic, flux_viscous,
                              mesh::P4estMesh{2},
-                             equations_parabolic::AbstractEquationsParabolic,
+                             equations_parabolic::AbstractEquationsParabolic{2},
                              surface_integral, dg::DG, cache)
     (; interfaces) = cache_parabolic
     (; contravariant_vectors) = cache_parabolic.elements
@@ -537,6 +851,98 @@ function prolong2interfaces!(cache_parabolic, flux_viscous,
     return nothing
 end
 
+# This is the version used when calculating the divergence of the viscous fluxes
+# of a 3D equation on a 2D p4est mesh
+# We pass the `surface_integral` argument solely for dispatch
+function prolong2interfaces!(cache_parabolic, flux_viscous,
+                             mesh::P4estMesh{2},
+                             equations_parabolic::AbstractEquationsParabolic{3},
+                             surface_integral, dg::DG, cache)
+    (; interfaces) = cache_parabolic
+    (; contravariant_vectors) = cache_parabolic.elements
+    index_range = eachnode(dg)
+    flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous
+
+    @threaded for interface in eachinterface(dg, cache)
+        # Copy solution data from the primary element using "delayed indexing" with
+        # a start value and a step size to get the correct face and orientation.
+        # Note that in the current implementation, the interface will be
+        # "aligned at the primary element", i.e., the index of the primary side
+        # will always run forwards.
+        primary_element = interfaces.neighbor_ids[1, interface]
+        primary_indices = interfaces.node_indices[1, interface]
+        primary_direction = indices2direction(primary_indices)
+
+        i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1],
+                                                                 index_range)
+        j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2],
+                                                                 index_range)
+
+        i_primary = i_primary_start
+        j_primary = j_primary_start
+        for i in eachnode(dg)
+
+            # this is the outward normal direction on the primary element
+            normal_direction = get_normal_direction(primary_direction,
+                                                    contravariant_vectors,
+                                                    i_primary, j_primary, primary_element)
+
+            for v in eachvariable(equations_parabolic)
+                # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*!
+                flux_viscous = SVector(flux_viscous_x[v, i_primary, j_primary,
+                                                      primary_element],
+                                       flux_viscous_y[v, i_primary, j_primary,
+                                                      primary_element],
+                                       flux_viscous_z[v, i_primary, j_primary,
+                                                      primary_element])
+
+                interfaces.u[1, v, i, interface] = dot(flux_viscous, normal_direction)
+            end
+            i_primary += i_primary_step
+            j_primary += j_primary_step
+        end
+
+        # Copy solution data from the secondary element using "delayed indexing" with
+        # a start value and a step size to get the correct face and orientation.
+        secondary_element = interfaces.neighbor_ids[2, interface]
+        secondary_indices = interfaces.node_indices[2, interface]
+        secondary_direction = indices2direction(secondary_indices)
+
+        i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1],
+                                                                     index_range)
+        j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2],
+                                                                     index_range)
+
+        i_secondary = i_secondary_start
+        j_secondary = j_secondary_start
+        for i in eachnode(dg)
+            # This is the outward normal direction on the secondary element.
+            # Here, we assume that normal_direction on the secondary element is
+            # the negative of normal_direction on the primary element.
+            normal_direction = get_normal_direction(secondary_direction,
+                                                    contravariant_vectors,
+                                                    i_secondary, j_secondary,
+                                                    secondary_element)
+
+            for v in eachvariable(equations_parabolic)
+                # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*!
+                flux_viscous = SVector(flux_viscous_x[v, i_secondary, j_secondary,
+                                                      secondary_element],
+                                       flux_viscous_y[v, i_secondary, j_secondary,
+                                                      secondary_element],
+                                       flux_viscous_z[v, i_secondary, j_secondary,
+                                                      secondary_element])
+                # store the normal flux with respect to the primary normal direction
+                interfaces.u[2, v, i, interface] = -dot(flux_viscous, normal_direction)
+            end
+            i_secondary += i_secondary_step
+            j_secondary += j_secondary_step
+        end
+    end
+
+    return nothing
+end
+
 # This version is used for divergence flux computations 
 function calc_interface_flux!(surface_flux_values,
                               mesh::P4estMesh{2}, equations_parabolic,
diff --git a/src/solvers/dgsem_tree/container_viscous_2d.jl b/src/solvers/dgsem_tree/container_viscous_2d.jl
index bd7ff413af5..436f6ae8c16 100644
--- a/src/solvers/dgsem_tree/container_viscous_2d.jl
+++ b/src/solvers/dgsem_tree/container_viscous_2d.jl
@@ -67,3 +67,77 @@ function Base.resize!(viscous_container::ViscousContainer2D, equations, dg, cach
     end
     return nothing
 end
+
+# Container for solving a 3D equation in a 2D grid
+mutable struct ViscousContainer2DEq3D{uEltype <: Real}
+    u_transformed::Array{uEltype, 4}
+    # Using an outer fixed-size datastructure leads to nasty implementations,
+    # see https://github.com/trixi-framework/Trixi.jl/pull/1629#discussion_r1355293953.
+    # Also: This does not result in speed up compared to using tuples for the internal 
+    # datastructures, see 
+    # https://github.com/trixi-framework/Trixi.jl/pull/1629#discussion_r1363352188.
+    gradients::Vector{Array{uEltype, 4}}
+    flux_viscous::Vector{Array{uEltype, 4}}
+
+    # internal `resize!`able storage
+    _u_transformed::Vector{uEltype}
+    # Use Tuple for outer, fixed-size datastructure
+    _gradients::Tuple{Vector{uEltype}, Vector{uEltype}, Vector{uEltype}}
+    _flux_viscous::Tuple{Vector{uEltype}, Vector{uEltype}, Vector{uEltype}}
+
+    function ViscousContainer2DEq3D{uEltype}(n_vars::Integer, n_nodes::Integer,
+                                             n_elements::Integer) where {uEltype <: Real}
+        new(Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements),
+            [Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements) for _ in 1:3],
+            [Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements) for _ in 1:3],
+            Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements),
+            (Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements),
+             Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements),
+             Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements)),
+            (Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements),
+             Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements),
+             Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements)))
+    end
+end
+
+function init_viscous_container_2d_eq_3D(n_vars::Integer, n_nodes::Integer,
+                                         n_elements::Integer,
+                                   ::Type{uEltype}) where {uEltype <: Real}
+    return ViscousContainer2DEq3D{uEltype}(n_vars, n_nodes, n_elements)
+end
+
+# Only one-dimensional `Array`s are `resize!`able in Julia.
+# Hence, we use `Vector`s as internal storage and `resize!`
+# them whenever needed. Then, we reuse the same memory by
+# `unsafe_wrap`ping multi-dimensional `Array`s around the
+# internal storage.
+# TODO: This resize can be merged with the one for ViscousContainer2D.
+function Base.resize!(viscous_container::ViscousContainer2DEq3D, equations, dg, cache)
+    capacity = nvariables(equations) * nnodes(dg) * nnodes(dg) * nelements(dg, cache)
+    resize!(viscous_container._u_transformed, capacity)
+    for dim in 1:3
+        resize!(viscous_container._gradients[dim], capacity)
+        resize!(viscous_container._flux_viscous[dim], capacity)
+    end
+
+    viscous_container.u_transformed = unsafe_wrap(Array,
+                                                  pointer(viscous_container._u_transformed),
+                                                  (nvariables(equations),
+                                                   nnodes(dg), nnodes(dg),
+                                                   nelements(dg, cache)))
+
+    for dim in 1:3
+        viscous_container.gradients[dim] = unsafe_wrap(Array,
+                                                       pointer(viscous_container._gradients[dim]),
+                                                       (nvariables(equations),
+                                                        nnodes(dg), nnodes(dg),
+                                                        nelements(dg, cache)))
+
+        viscous_container.flux_viscous[dim] = unsafe_wrap(Array,
+                                                          pointer(viscous_container._flux_viscous[dim]),
+                                                          (nvariables(equations),
+                                                           nnodes(dg), nnodes(dg),
+                                                           nelements(dg, cache)))
+    end
+    return nothing
+end
diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
index 3083ae30680..5e4186678d4 100644
--- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
@@ -294,7 +294,7 @@ end
 function calc_viscous_fluxes!(flux_viscous,
                               gradients, u_transformed,
                               mesh::Union{TreeMesh{2}, P4estMesh{2}},
-                              equations_parabolic::AbstractEquationsParabolic,
+                              equations_parabolic::AbstractEquationsParabolic{2},
                               dg::DG, cache, cache_parabolic)
     gradients_x, gradients_y = gradients
     flux_viscous_x, flux_viscous_y = flux_viscous # output arrays
@@ -322,6 +322,43 @@ function calc_viscous_fluxes!(flux_viscous,
     end
 end
 
+function calc_viscous_fluxes!(flux_viscous,
+                              gradients, u_transformed,
+                              mesh::Union{TreeMesh{2}, P4estMesh{2}},
+                              equations_parabolic::AbstractEquationsParabolic{3},
+                              dg::DG, cache, cache_parabolic)
+    gradients_x, gradients_y, gradients_z = gradients
+    flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous # output arrays
+
+    @threaded for element in eachelement(dg, cache)
+        for j in eachnode(dg), i in eachnode(dg)
+            # Get solution and gradients
+            u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, j,
+                                   element)
+            gradients_1_node = get_node_vars(gradients_x, equations_parabolic, dg, i, j,
+                                             element)
+            gradients_2_node = get_node_vars(gradients_y, equations_parabolic, dg, i, j,
+                                             element)
+            gradients_3_node = get_node_vars(gradients_z, equations_parabolic, dg, i, j,
+                                             element)
+            
+            # Calculate viscous flux and store each component for later use
+            flux_viscous_node_x = flux(u_node, (gradients_1_node, gradients_2_node, gradients_3_node), 1,
+                                       equations_parabolic)
+            flux_viscous_node_y = flux(u_node, (gradients_1_node, gradients_2_node, gradients_3_node), 2,
+                                       equations_parabolic)
+            flux_viscous_node_z = flux(u_node, (gradients_1_node, gradients_2_node, gradients_3_node), 3,
+                                       equations_parabolic)
+            set_node_vars!(flux_viscous_x, flux_viscous_node_x, equations_parabolic, dg,
+                           i, j, element)
+            set_node_vars!(flux_viscous_y, flux_viscous_node_y, equations_parabolic, dg,
+                           i, j, element)
+            set_node_vars!(flux_viscous_z, flux_viscous_node_z, equations_parabolic, dg,
+                           i, j, element)
+        end
+    end
+end
+
 # TODO: parabolic; decide if we should keep this, and if so, extend to 3D.
 function get_unsigned_normal_vector_2d(direction)
     if direction > 4 || direction < 1