diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index e85f783f1de6..8d935dc41441 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -4859,6 +4859,15 @@ class CConfig {
    */
   bool GetMUSCL_Flow(void) const { return MUSCL_Flow; }
 
+  /*! \brief Whether the turbulence equations reconstruct. */
+  bool GetMUSCL_Turb(void) const { return MUSCL_Turb; }
+
+  /*! \brief Whether the species equations reconstruct. */
+  bool GetMUSCL_Species(void) const { return MUSCL_Species; }
+
+  /*! \brief Whether the heat equation reconstructs. */
+  bool GetMUSCL_Heat(void) const { return MUSCL_Heat; }
+
   /*!
    * \brief Get if the upwind scheme used MUSCL or not.
    * \return MUSCL scheme.
diff --git a/SU2_CFD/include/numerics/scalar/scalar_edge_flux.hpp b/SU2_CFD/include/numerics/scalar/scalar_edge_flux.hpp
index b977b010b746..df1275c0ceeb 100644
--- a/SU2_CFD/include/numerics/scalar/scalar_edge_flux.hpp
+++ b/SU2_CFD/include/numerics/scalar/scalar_edge_flux.hpp
@@ -62,6 +62,12 @@ struct ScalarFluxOptions {
   bool viscous = false;   /*!< \brief Whether the diffusion term contributes. */
   bool oneSided = false;  /*!< \brief Whether only the row of i is assembled. */
   bool muscl = false;     /*!< \brief Whether the convective scheme reconstructs. */
+  /*!
+   * \brief Bounds of the transported variables (those that clip the solution). A reconstructed face
+   *        value outside them falls back to the value of its point (first order). Null: no bounds.
+   */
+  const su2double* lowerLimit = nullptr;
+  const su2double* upperLimit = nullptr;
 
   /*!
    * \brief Options of the interior edge loop: both terms, both rows, and reconstruction and
@@ -349,10 +355,9 @@ class CUpwScalarBase : public CUpwScalarFlux<Double_, Derived, FlowIndices, nDim
 
   /*!
    * \brief MUSCL reconstruction parameters, read from CConfig once per construction (i.e. once
-   *        per nonlinear iteration, see CScalarSolver::EdgeFluxResidual) instead of per edge;
-   *        this is also where the scalar limiter's freezing (GetLimiterIter) is resolved, by
-   *        collapsing its type to NONE once frozen. The flow limiter is not frozen this way: once
-   *        the flow solver stops recomputing it, it keeps applying the last values it has.
+   *        per nonlinear iteration, see CScalarSolver::EdgeFluxResidual) instead of per edge.
+   *        Point limiters keep applying their last stored values after LIMITER_ITER.
+   *        VAN_ALBADA_EDGE is computed on every edge and is never frozen.
    */
   const su2double kappa, umusclRamp, kappaFlow;
   const LIMITER limiterType, limiterTypeFlow;
@@ -372,9 +377,8 @@ class CUpwScalarBase : public CUpwScalarFlux<Double_, Derived, FlowIndices, nDim
         kappa(config.GetMUSCL_Kappa()),
         umusclRamp(config.GetMUSCLRampValue()),
         kappaFlow(config.GetMUSCL_Kappa_Flow()),
-        limiterType(config.GetInnerIter() <= config.GetLimiterIter() ? config.GetKind_SlopeLimit() : LIMITER::NONE),
-        limiterTypeFlow(config.GetKind_SlopeLimit_Flow() != LIMITER::VAN_ALBADA_EDGE ? config.GetKind_SlopeLimit_Flow()
-                                                                                     : LIMITER::NONE),
+        limiterType(config.GetKind_SlopeLimit()),
+        limiterTypeFlow(config.GetKind_SlopeLimit_Flow()),
         musclFlow(config.GetMUSCL_Flow() && config.GetKind_ConvNumScheme_Flow() == SPACE_UPWIND) {
     if (nEqn > Size) {
       SU2_MPI::Error("The number of equations exceeds MAX_TRANSPORTED_SPECIES, the size of the static arrays. Increase it\n"
@@ -474,14 +478,29 @@ class CUpwScalarBase : public CUpwScalarFlux<Double_, Derived, FlowIndices, nDim
       }
 
       if (opt.muscl) {
+        const auto phi1st = phi;
         if constexpr (nVar != Dynamic) {
-          reconstruct<nVar>(iPoint, jPoint, vector_ij, side_i.scalarNodes.GetGradient_Reconstruction(),
-                            side_i.scalarNodes.GetLimiter(), limiterType, 0, phi, kappa, umusclRamp);
+          reconstruct<nVar, MusclLimiter::EDGE_SCALAR>(
+              iPoint, jPoint, vector_ij, side_i.scalarNodes.GetGradient_Reconstruction(),
+              side_i.scalarNodes.GetLimiter(), limiterType, 0, phi, kappa, umusclRamp);
         } else {
           /*--- A dynamic model's equation count is only known at runtime, so the reconstructed
            * width is passed as an argument instead of a template parameter. ---*/
-          reconstruct(iPoint, jPoint, vector_ij, side_i.scalarNodes.GetGradient_Reconstruction(),
-                      side_i.scalarNodes.GetLimiter(), limiterType, 0, phi, kappa, umusclRamp, res.nVar);
+          reconstruct<0, MusclLimiter::EDGE_SCALAR>(
+              iPoint, jPoint, vector_ij, side_i.scalarNodes.GetGradient_Reconstruction(),
+              side_i.scalarNodes.GetLimiter(), limiterType, 0, phi, kappa, umusclRamp, res.nVar);
+        }
+
+        /*--- A face value outside the bounds of the variable (e.g. a negative k or omega next to a wall,
+         * where omega varies by orders of magnitude) falls back to the value of its point. ---*/
+        if (opt.lowerLimit) {
+          for (size_t iVar = 0; iVar < res.nVar; ++iVar) {
+            const su2double lo = opt.lowerLimit[iVar], hi = opt.upperLimit[iVar];
+            const Double out_i = fmax(phi.i.all(iVar) < lo, phi.i.all(iVar) > hi);
+            const Double out_j = fmax(phi.j.all(iVar) < lo, phi.j.all(iVar) > hi);
+            phi.i.all(iVar) = out_i * phi1st.i.all(iVar) + (1 - out_i) * phi.i.all(iVar);
+            phi.j.all(iVar) = out_j * phi1st.j.all(iVar) + (1 - out_j) * phi.j.all(iVar);
+          }
         }
       }
 
diff --git a/SU2_CFD/include/numerics/util.hpp b/SU2_CFD/include/numerics/util.hpp
index cb617e5dd22c..937e403d197e 100644
--- a/SU2_CFD/include/numerics/util.hpp
+++ b/SU2_CFD/include/numerics/util.hpp
@@ -34,6 +34,8 @@
 #include "../../../Common/include/containers/C2DContainer.hpp"
 #include "../../../Common/include/linear_algebra/CSysVector.hpp"
 #include "../../../Common/include/linear_algebra/CSysMatrix.hpp"
+#include "../../../Common/include/geometry/CGeometry.hpp"
+#include "../limiters/CLimiterDetails.hpp"
 
 /*!
  * \enum UpdateType
@@ -498,8 +500,11 @@ FORCEINLINE auto gradientReader(typename CLaneTraits<Double>::Int iPoint, const
 
 /*!
  * \brief How the reconstructed differences are limited.
+ * \note EDGE is the van Albada edge limiter of the flow solver, which the scalar solvers also use to
+ *       reconstruct the flow velocity so that both see the same face velocity. EDGE_SCALAR is the
+ *       van Albada function of the point limiters, used for the transported scalar variables.
  */
-enum class MusclLimiter { NONE, EDGE, POINT };
+enum class MusclLimiter { NONE, EDGE, EDGE_SCALAR, POINT };
 
 /*!
  * \brief U-MUSCL reconstruction of nVarGrad variables, from the gradient rows starting at iRow.
@@ -534,6 +539,9 @@ FORCEINLINE void muscl(typename CLaneTraits<Double>::Int iPoint, typename CLaneT
       const Double delta_ij_2 = pow(delta_ij, 2) + 1e-6;
       lim_i = (delta_ij_2 + proj_i * delta_ij) / (pow(proj_i, 2) + delta_ij_2);
       lim_j = (delta_ij_2 + proj_j * delta_ij) / (pow(proj_j, 2) + delta_ij_2);
+    } else if constexpr (limiterKind == MusclLimiter::EDGE_SCALAR) {
+      lim_i = LimiterHelpers<Double>::vanAlbadaFunction(proj_i, delta_ij, Double(EPS));
+      lim_j = LimiterHelpers<Double>::vanAlbadaFunction(proj_j, delta_ij, Double(EPS));
     } else if constexpr (limiterKind == MusclLimiter::POINT) {
       lim_i = gatherVariables(iPoint, limiter, iRow + iVar);
       lim_j = gatherVariables(jPoint, limiter, iRow + iVar);
@@ -548,21 +556,25 @@ FORCEINLINE void muscl(typename CLaneTraits<Double>::Int iPoint, typename CLaneT
 /*!
  * \brief Reconstruct a slice of nVarGrad variables starting at column iRow, dispatching on the
  *        limiter type; shared by the flow and the scalar reconstructions so both call one body.
+ * \note edgeLimiter is the form that VAN_ALBADA_EDGE uses (EDGE for flow variables, EDGE_SCALAR for
+ *       transported scalars).
  */
-template <size_t nVarGrad_ = 0, size_t nDim, class Double, class VarType, class Limiter_t, class Gradient_t>
+template <size_t nVarGrad_ = 0, MusclLimiter edgeLimiter = MusclLimiter::EDGE, size_t nDim, class Double, class VarType,
+          class Limiter_t, class Gradient_t>
 FORCEINLINE void reconstruct(typename CLaneTraits<Double>::Int iPoint, typename CLaneTraits<Double>::Int jPoint,
                              const Vector<Double, nDim>& vector_ij, const Gradient_t& gradient,
                              const Limiter_t& limiter, LIMITER limiterType, size_t iRow, CPair<VarType>& V,
                              const CNonDeduced<Double>& kappa, const CNonDeduced<Double>& umusclRamp,
                              size_t nVarGradRuntime = 0) {
+  static_assert(edgeLimiter == MusclLimiter::EDGE || edgeLimiter == MusclLimiter::EDGE_SCALAR);
   switch (limiterType) {
     case LIMITER::NONE:
       muscl<MusclLimiter::NONE, nVarGrad_>(iPoint, jPoint, vector_ij, gradient, limiter, iRow, V, kappa, umusclRamp,
                                            nVarGradRuntime);
       break;
     case LIMITER::VAN_ALBADA_EDGE:
-      muscl<MusclLimiter::EDGE, nVarGrad_>(iPoint, jPoint, vector_ij, gradient, limiter, iRow, V, kappa, umusclRamp,
-                                           nVarGradRuntime);
+      muscl<edgeLimiter, nVarGrad_>(iPoint, jPoint, vector_ij, gradient, limiter, iRow, V, kappa, umusclRamp,
+                                    nVarGradRuntime);
       break;
     default:
       muscl<MusclLimiter::POINT, nVarGrad_>(iPoint, jPoint, vector_ij, gradient, limiter, iRow, V, kappa, umusclRamp,
diff --git a/SU2_CFD/include/solvers/CScalarSolver.hpp b/SU2_CFD/include/solvers/CScalarSolver.hpp
index 50ecf8a42c7f..a7a559527fbc 100644
--- a/SU2_CFD/include/solvers/CScalarSolver.hpp
+++ b/SU2_CFD/include/solvers/CScalarSolver.hpp
@@ -190,7 +190,8 @@ class CScalarSolver : public CSolver {
           if (conv_numerics->GetBoundedScalar()) {
             const su2double* velocity = &PrimVar_j[prim_idx.Velocity()];
             const su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint);
-            conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, true, density, velocity, Normal));
+            conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, true, density, velocity, Normal,
+                                                         &PrimVar_j[prim_idx.Density()], weight));
           }
 
           auto residual = conv_numerics->ComputeResidual(config);
@@ -244,16 +245,35 @@ class CScalarSolver : public CSolver {
    * \brief Applies a convective flux correction to negate the effects of flow divergence at a BC node.
    * \note This function should be used for nodes that are part of a boundary marker, it computes a mass flux
    * from density and velocity at the node, and the outward-pointing normal (-1 * normal of vertex).
+   * \param[in] weight - Weight of the boundary state, e.g. of a donor at a sliding interface, which also weights
+   *            the flux (the returned mass flux is not weighted).
    * \return The mass flux.
    */
   inline su2double BoundedScalarBCFlux(unsigned long iPoint, bool implicit, const su2double& density,
-                                       const su2double* velocity, const su2double* normal) {
-    const su2double edgeMassFlux = density * GeometryToolbox::DotProduct(nDim, velocity, normal);
-    LinSysRes.AddBlock(iPoint, nodes->GetSolution(iPoint), -edgeMassFlux);
-    if (implicit) Jacobian.AddVal2Diag(iPoint, -edgeMassFlux);
+                                       const su2double* velocity, const su2double* normal,
+                                       const su2double* densityGhost = nullptr, su2double weight = 1.0) {
+    const su2double rhoGhost = densityGhost ? *densityGhost : density;
+    const su2double volumeFlux = GeometryToolbox::DotProduct(nDim, velocity, normal);
+    const su2double edgeMassFlux = volumeFlux * (volumeFlux > 0.0 ? density : rhoGhost);
+    const su2double q = weight * BoundedScalarDivergenceFlux(edgeMassFlux, density, rhoGhost);
+    LinSysRes.AddBlock(iPoint, nodes->GetSolution(iPoint), -q);
+    if (implicit) Jacobian.AddVal2Diag(iPoint, -q / (Conservative ? density : 1.0));
     return edgeMassFlux;
   }
 
+  /*!
+   * \brief Face flux of the bounded scalar divergence correction, -q * phi_i, such that it cancels the
+   *        part of the upwind flux due to the discrete divergence of the mass flux. It uses the same
+   *        coefficients as the flux: the mass flux for a conservative scalar (whose flux is
+   *        mdot+ phi_i + mdot- phi_j), the upwinded volume flux mdot+/rho_i + mdot-/rho_j otherwise.
+   * \note The Jacobian w.r.t. the conserved variable of a conservative scalar, rho * phi, is -q / rho_i.
+   */
+  inline su2double BoundedScalarDivergenceFlux(const su2double& massFlux, const su2double& rho_i,
+                                               const su2double& rho_j) const {
+    if (Conservative) return massFlux;
+    return fmax(massFlux, 0.0) / rho_i + fmin(massFlux, 0.0) / rho_j;
+  }
+
   /*!
    * \brief Move solution to previous time levels (for restarts).
    */
diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl
index af03659fe7f5..8d68faf26a77 100644
--- a/SU2_CFD/include/solvers/CScalarSolver.inl
+++ b/SU2_CFD/include/solvers/CScalarSolver.inl
@@ -118,7 +118,9 @@ void CScalarSolver<VariableType>::CommonPreprocessing(CGeometry *geometry, const
    * before calling these solver functions. ---*/
   const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
   const bool muscl = config->GetMUSCL();
+  /*--- The edge limiter (VAN_ALBADA_EDGE) is applied in the flux kernel and needs no point values. ---*/
   const bool limiter = (config->GetKind_SlopeLimit() != LIMITER::NONE) &&
+                       (config->GetKind_SlopeLimit() != LIMITER::VAN_ALBADA_EDGE) &&
                        (config->GetInnerIter() <= config->GetLimiterIter());
 
   /*--- Clear residual and system matrix, not needed for
@@ -132,9 +134,11 @@ void CScalarSolver<VariableType>::CommonPreprocessing(CGeometry *geometry, const
     }
   }
 
-  /*--- Upwind second order reconstruction and gradients ---*/
+  /*--- Upwind second order reconstruction and gradients. The reconstruction gradient is global to the
+   * case (it is needed as soon as any equation system reconstructs), so it is skipped here when this
+   * system does not reconstruct: only the reconstruction and the limiter read it. ---*/
 
-  if (config->GetReconstructionGradientRequired()) {
+  if (muscl && config->GetReconstructionGradientRequired()) {
     switch(config->GetKind_Gradient_Method_Recon()) {
       case GREEN_GAUSS: SetSolution_Gradient_GG(geometry, config, -1, true); break;
       case LEAST_SQUARES: SetSolution_Gradient_LS(geometry, config, -1, true); break;
@@ -174,9 +178,14 @@ void CScalarSolver<VariableType>::SumEdgeFluxes(const CGeometry* geometry) {
 template <class VariableType>
 template <class Scheme>
 void CScalarSolver<VariableType>::EdgeFluxResidual(const CGeometry* geometry, CSolver** solver_container,
-                                                    const CConfig* config, const ScalarFluxOptions& opt) {
+                                                    const CConfig* config, const ScalarFluxOptions& optIn) {
   SU2_ZONE_SCOPED
 
+  /*--- The reconstructed face values are kept within the bounds that clip the solution. ---*/
+  auto opt = optIn;
+  opt.lowerLimit = lowerlimit;
+  opt.upperLimit = upperlimit;
+
   using Double = typename Scheme::Double;
   constexpr int nDim = Scheme::nDim;
 
@@ -217,11 +226,13 @@ void CScalarSolver<VariableType>::EdgeFluxResidual(const CGeometry* geometry, CS
       /*--- Bounded scalar divergence correction, per edge; the ReducerStrategy equivalent runs
        * in a per-point pass below, where the diagonal is not written from the edge loop. ---*/
       if (opt.boundedScalar && !ReducerStrategy) {
-        LinSysRes.AddBlock(iPoint, nodes->GetSolution(iPoint), -massFlux);
-        LinSysRes.AddBlock(jPoint, nodes->GetSolution(jPoint), massFlux);
+        const su2double rho_i = flowNodes->GetDensity(iPoint), rho_j = flowNodes->GetDensity(jPoint);
+        const su2double q = BoundedScalarDivergenceFlux(massFlux, rho_i, rho_j);
+        LinSysRes.AddBlock(iPoint, nodes->GetSolution(iPoint), -q);
+        LinSysRes.AddBlock(jPoint, nodes->GetSolution(jPoint), q);
         if (opt.implicit) {
-          Jacobian.AddVal2Diag(iPoint, -massFlux);
-          Jacobian.AddVal2Diag(jPoint, massFlux);
+          Jacobian.AddVal2Diag(iPoint, -q / (Conservative ? rho_i : 1.0));
+          Jacobian.AddVal2Diag(jPoint, q / (Conservative ? rho_j : 1.0));
         }
       }
     }
@@ -239,15 +250,19 @@ void CScalarSolver<VariableType>::EdgeFluxResidual(const CGeometry* geometry, CS
       SU2_OMP_FOR_STAT(omp_chunk_size)
       for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) {
         const auto* solution = nodes->GetSolution(iPoint);
+        const su2double rho_i = flowNodes->GetDensity(iPoint);
         su2double divergence = 0;
 
-        for (auto iEdge : geometry->nodes->GetEdges(iPoint)) {
+        for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) {
+          const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh);
+          const auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh);
           const auto sign = (iPoint == geometry->edges->GetNode(iEdge, 0)) ? 1 : -1;
-          const su2double edgeMassFlux = sign * (*edgeMassFluxes)[iEdge];
-          divergence += edgeMassFlux;
-          LinSysRes.AddBlock(iPoint, solution, -edgeMassFlux);
+          const su2double q = BoundedScalarDivergenceFlux(sign * (*edgeMassFluxes)[iEdge], rho_i,
+                                                          flowNodes->GetDensity(jPoint));
+          divergence += q;
+          LinSysRes.AddBlock(iPoint, solution, -q);
         }
-        if (opt.implicit) Jacobian.AddVal2Diag(iPoint, -divergence);
+        if (opt.implicit) Jacobian.AddVal2Diag(iPoint, -divergence / (Conservative ? rho_i : 1.0));
       }
       END_SU2_OMP_FOR
     }
@@ -282,8 +297,9 @@ void CScalarSolver<VariableType>::BoundaryFluxResidual(const CGeometry* geometry
 
     Double massFlux = 0.0;
     if (opt.boundedScalar) {
+      const auto* ghostPrim = ghostFlowNodes->GetPrimitive(iVertex);
       massFlux = BoundedScalarBCFlux(iPoint, opt.implicit, flowNodes->GetDensity(iPoint),
-                                     &ghostFlowNodes->GetPrimitive(iVertex)[prim_idx.Velocity()], normal.data());
+                                     &ghostPrim[prim_idx.Velocity()], normal.data(), &ghostPrim[prim_idx.Density()]);
     }
 
     const auto res = flux.ComputeFlux(opt, iPoint, side_i, iVertex, side_j, normal, massFlux);
@@ -347,7 +363,8 @@ void CScalarSolver<VariableType>::FluidInterfaceFluxResidual(const CGeometry* ge
         su2double massFlux = 0.0;
         if (optConv.boundedScalar) {
           massFlux = BoundedScalarBCFlux(iPoint, optConv.implicit, flowNodes->GetDensity(iPoint),
-                                         &PrimVar_j[prim_idx.Velocity()], normal.data());
+                                         &PrimVar_j[prim_idx.Velocity()], normal.data(), &PrimVar_j[prim_idx.Density()],
+                                         weight);
         }
 
         const auto res = flux.ComputeFlux(optConv, iPoint, side_i, iVertex, side_j, normal, massFlux);
@@ -446,7 +463,8 @@ void CScalarSolver<VariableType>::BC_Far_Field(CGeometry* geometry, CSolver** so
       if (conv_numerics->GetBoundedScalar()) {
         const su2double* velocity = &V_infty[prim_idx.Velocity()];
         const su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint);
-        conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, implicit, density, velocity, Normal));
+        conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, implicit, density, velocity, Normal,
+                                                     &V_infty[prim_idx.Density()]));
       }
 
       /*--- Compute residuals and Jacobians ---*/
diff --git a/SU2_CFD/include/variables/CScalarVariable.hpp b/SU2_CFD/include/variables/CScalarVariable.hpp
index 9c9540673ad8..380e544312eb 100644
--- a/SU2_CFD/include/variables/CScalarVariable.hpp
+++ b/SU2_CFD/include/variables/CScalarVariable.hpp
@@ -49,8 +49,9 @@ class CScalarVariable : public CVariable {
    * \param[in] ndim - Number of dimensions of the problem.
    * \param[in] nvar - Number of variables of the problem.
    * \param[in] config - Definition of the particular problem.
+   * \param[in] muscl - Whether this scalar system reconstructs.
    */
-  CScalarVariable(unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config);
+  CScalarVariable(unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config, bool muscl);
 
   /*!
    * \brief Get the array of the reconstruction variables gradient at a node.
diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp
index f454556c7ecd..189b9dc53104 100644
--- a/SU2_CFD/src/output/CFlowOutput.cpp
+++ b/SU2_CFD/src/output/CFlowOutput.cpp
@@ -1477,7 +1477,7 @@ void CFlowOutput::SetVolumeOutputFieldsScalarLimiter(const CConfig* config) {
   /*--- Only place outputs of the "SOLUTION" group for species transport here. ---*/
 
 
-  if (config->GetKind_SlopeLimit_Turb() != LIMITER::NONE) {
+  if ((config->GetKind_SlopeLimit_Turb() != LIMITER::NONE && config->GetKind_SlopeLimit_Turb() != LIMITER::VAN_ALBADA_EDGE)) {
     switch (TurbModelFamily(config->GetKind_Turb_Model())) {
       case TURB_FAMILY::SA:
         AddVolumeOutput("LIMITER_NU_TILDE", "Limiter_Nu_Tilde", "LIMITER", "Limiter value of the Spalart-Allmaras variable");
@@ -1493,7 +1493,7 @@ void CFlowOutput::SetVolumeOutputFieldsScalarLimiter(const CConfig* config) {
     }
   }
 
-  if (config->GetKind_SlopeLimit_Species() != LIMITER::NONE) {
+  if ((config->GetKind_SlopeLimit_Species() != LIMITER::NONE && config->GetKind_SlopeLimit_Species() != LIMITER::VAN_ALBADA_EDGE)) {
     switch (config->GetKind_Species_Model()) {
       case SPECIES_MODEL::SPECIES_TRANSPORT:
         for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++)
@@ -1667,7 +1667,7 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
     SetVolumeOutputValue("Q_CRITERION", iPoint, GetQCriterion(Node_Flow->GetVelocityGradient(iPoint)));
   }
 
-  const bool limiter = (config->GetKind_SlopeLimit_Turb() != LIMITER::NONE);
+  const bool limiter = ((config->GetKind_SlopeLimit_Turb() != LIMITER::NONE && config->GetKind_SlopeLimit_Turb() != LIMITER::VAN_ALBADA_EDGE));
 
   switch (TurbModelFamily(config->GetKind_Turb_Model())) {
     case TURB_FAMILY::SA:
@@ -1748,10 +1748,11 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
         SetVolumeOutputValue("DIFFUSIVITY_" + label, iPoint, Node_Species->GetDiffusivity(iPoint, iVar));
         if (config->GetCanteraCombustion())
           SetVolumeOutputValue("CHEMICAL_SOURCE_TERM_" + label, iPoint, Node_Species->GetChemicalSourceTerm(iPoint, iVar));
-        if (config->GetKind_SlopeLimit_Species() != LIMITER::NONE)
+        if (config->GetKind_SlopeLimit_Species() != LIMITER::NONE && config->GetKind_SlopeLimit_Species() != LIMITER::VAN_ALBADA_EDGE)
           SetVolumeOutputValue("LIMITER_SPECIES_" + label, iPoint, Node_Species->GetLimiter(iPoint, iVar));
         if (config->GetPyCustomSource())
           SetVolumeOutputValue("SPECIES_UDS_" + label, iPoint, Node_Species->GetUserDefinedSource()(iPoint, iVar));
+
       }
       if (config->GetCanteraCombustion())
         SetVolumeOutputValue("HEAT_RELEASE", iPoint, Node_Species->GetHeatRelease(iPoint));
@@ -1778,7 +1779,7 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
         SetVolumeOutputValue("RES_" + scalar_name, iPoint, solver[SPECIES_SOL]->LinSysRes(iPoint, flamelet_config_options.n_control_vars + i_scalar));
       }
 
-      if (config->GetKind_SlopeLimit_Species() != LIMITER::NONE) {
+      if ((config->GetKind_SlopeLimit_Species() != LIMITER::NONE && config->GetKind_SlopeLimit_Species() != LIMITER::VAN_ALBADA_EDGE)) {
         /*--- Limiter for controlling variable transport equations. ---*/
         for (auto iCV=0u; iCV<flamelet_config_options.n_control_vars; iCV++) {
           const auto& cv_name = flamelet_config_options.controlling_variable_names[iCV];
diff --git a/SU2_CFD/src/variables/CHeatVariable.cpp b/SU2_CFD/src/variables/CHeatVariable.cpp
index 9bc30e42a15f..0060acc3df7d 100644
--- a/SU2_CFD/src/variables/CHeatVariable.cpp
+++ b/SU2_CFD/src/variables/CHeatVariable.cpp
@@ -28,7 +28,7 @@
 #include "../../include/variables/CHeatVariable.hpp"
 
 CHeatVariable::CHeatVariable(su2double heat, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config)
-  : CScalarVariable(npoint, ndim, nvar, config) {
+  : CScalarVariable(npoint, ndim, nvar, config, config->GetMUSCL_Heat()) {
 
   const bool dual_time = ((config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) ||
                           (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND));
diff --git a/SU2_CFD/src/variables/CPoissonVariable.cpp b/SU2_CFD/src/variables/CPoissonVariable.cpp
index 0ac42afb417b..0287afb0d0e3 100644
--- a/SU2_CFD/src/variables/CPoissonVariable.cpp
+++ b/SU2_CFD/src/variables/CPoissonVariable.cpp
@@ -28,7 +28,7 @@
 #include "../../include/variables/CPoissonVariable.hpp"
 
 CPoissonVariable::CPoissonVariable(su2double value, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config)
-  : CScalarVariable(npoint, ndim, nvar, config) {
+  : CScalarVariable(npoint, ndim, nvar, config, false) {
 
   /*--- Initialization ---*/
 
diff --git a/SU2_CFD/src/variables/CScalarVariable.cpp b/SU2_CFD/src/variables/CScalarVariable.cpp
index b7437011cd24..9ea4eb7acd08 100644
--- a/SU2_CFD/src/variables/CScalarVariable.cpp
+++ b/SU2_CFD/src/variables/CScalarVariable.cpp
@@ -27,14 +27,15 @@
 
 #include "../../include/variables/CScalarVariable.hpp"
 
-CScalarVariable::CScalarVariable(unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config)
+CScalarVariable::CScalarVariable(unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config,
+                                 bool muscl)
     : CVariable(npoint, ndim, nvar, config),
-      Gradient_Reconstruction(config->GetReconstructionGradientRequired() ? Gradient_Aux : Gradient) {
+      Gradient_Reconstruction(muscl && config->GetReconstructionGradientRequired() ? Gradient_Aux : Gradient) {
   /*--- Gradient related fields ---*/
 
   Gradient.resize(nPoint, nVar, nDim, 0.0);
 
-  if (config->GetReconstructionGradientRequired()) {
+  if (muscl && config->GetReconstructionGradientRequired()) {
     Gradient_Aux.resize(nPoint, nVar, nDim, 0.0);
   }
 
diff --git a/SU2_CFD/src/variables/CSpeciesVariable.cpp b/SU2_CFD/src/variables/CSpeciesVariable.cpp
index 103ff2505c27..ab46054411e8 100644
--- a/SU2_CFD/src/variables/CSpeciesVariable.cpp
+++ b/SU2_CFD/src/variables/CSpeciesVariable.cpp
@@ -29,8 +29,9 @@
 
 CSpeciesVariable::CSpeciesVariable(const su2double* species_inf, unsigned long npoint, unsigned long ndim,
                                    unsigned long nvar, const CConfig* config)
-    : CScalarVariable(npoint, ndim, nvar, config) {
+    : CScalarVariable(npoint, ndim, nvar, config, config->GetMUSCL_Species()) {
   /*--- Allocate space for the mass diffusivity and chemical source term. ---*/
+
   Diffusivity.resize(nPoint, nVar + 1) = su2double(0.0);
   if (config->GetCanteraCombustion()) {
     SpeciesSourceTerm.resize(nPoint, nVar + 1) = su2double(0.0);
diff --git a/SU2_CFD/src/variables/CTurbVariable.cpp b/SU2_CFD/src/variables/CTurbVariable.cpp
index 139223055a42..71ca449ba920 100644
--- a/SU2_CFD/src/variables/CTurbVariable.cpp
+++ b/SU2_CFD/src/variables/CTurbVariable.cpp
@@ -30,7 +30,7 @@
 
 
 CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config)
-  : CScalarVariable(npoint, ndim, nvar, config) {
+  : CScalarVariable(npoint, ndim, nvar, config, config->GetMUSCL_Turb()) {
 
     turb_index.resize(nPoint) = su2double(1.0);
     intermittency.resize(nPoint) = su2double(1.0);
diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py
index 21c63694f567..d6149bafbbba 100644
--- a/TestCases/hybrid_regression.py
+++ b/TestCases/hybrid_regression.py
@@ -170,7 +170,7 @@ def main():
     rae2822_sa.cfg_dir   = "rans/rae2822"
     rae2822_sa.cfg_file  = "turb_SA_RAE2822.cfg"
     rae2822_sa.test_iter = 20
-    rae2822_sa.test_vals = [-2.916839, -5.161263, 0.802145, 0.019586, 1000.000000]
+    rae2822_sa.test_vals = [-2.915770, -5.310335, 0.801323, 0.019588, 1000.000000]
     test_list.append(rae2822_sa)
 
     # RAE2822 SST
@@ -178,7 +178,7 @@ def main():
     rae2822_sst.cfg_dir   = "rans/rae2822"
     rae2822_sst.cfg_file  = "turb_SST_RAE2822.cfg"
     rae2822_sst.test_iter = 20
-    rae2822_sst.test_vals = [-1.461774, 5.883885, 0.630267, 0.020404, 100.000000]
+    rae2822_sst.test_vals = [-1.472874, 5.886555, 0.669542, 0.023434, 100.000000]
     test_list.append(rae2822_sst)
 
     # RAE2822 SST_SUST
@@ -452,7 +452,7 @@ def main():
     inc_turb_naca0012.cfg_dir   = "incomp_rans/naca0012"
     inc_turb_naca0012.cfg_file  = "naca0012.cfg"
     inc_turb_naca0012.test_iter = 20
-    inc_turb_naca0012.test_vals = [-4.758114, -10.974548, -0.000006, -0.028637, 4.000000, -4.016988, 2.000000, -4.487046]
+    inc_turb_naca0012.test_vals = [-4.758114, -10.974548, -0.000006, -0.028637, 4.000000, -4.016947, 2.000000, -4.490161]
     test_list.append(inc_turb_naca0012)
 
     # NACA0012, SST_SUST
@@ -711,7 +711,7 @@ def main():
     bars_SST_2D.cfg_dir   = "sliding_interface/bars_SST_2D"
     bars_SST_2D.cfg_file  = "bars.cfg"
     bars_SST_2D.test_iter = 13
-    bars_SST_2D.test_vals = [13.000000, -0.455981, -1.541069]
+    bars_SST_2D.test_vals = [13.000000, -0.462453, -1.541069]
     bars_SST_2D.multizone = True
     test_list.append(bars_SST_2D)
 
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index 89bdeb0fe118..013464c7fe76 100755
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -78,7 +78,8 @@ def main():
     cfd_flamelet_ch4_axi.cfg_dir = "flamelet/05_laminar_premixed_ch4_flame_cfd_axi"
     cfd_flamelet_ch4_axi.cfg_file = "lam_prem_ch4_cfd_axi.cfg"
     cfd_flamelet_ch4_axi.test_iter = 10
-    cfd_flamelet_ch4_axi.test_vals = [-9.774148, -11.255837, -10.017962, -11.850828, -4.916261, 72.669000]
+    cfd_flamelet_ch4_axi.test_vals = [-9.774148, -11.256092, -10.017960, -11.850827, -4.916267, 72.669000]
+
     cfd_flamelet_ch4_axi.new_output = True
     test_list.append(cfd_flamelet_ch4_axi)
 
@@ -113,7 +114,7 @@ def main():
     cfd_flamelet_ch4_unsteady.cfg_dir = "flamelet/09_laminar_premixed_ch4_flame_unsteady"
     cfd_flamelet_ch4_unsteady.cfg_file = "lam_prem_ch4_unsteady.cfg"
     cfd_flamelet_ch4_unsteady.test_iter = 5
-    cfd_flamelet_ch4_unsteady.test_vals = [-8.856417, -8.095249, -9.153745, -9.321678]
+    cfd_flamelet_ch4_unsteady.test_vals = [-8.856677, -8.095263, -9.153743, -9.321676]
     cfd_flamelet_ch4_unsteady.test_vals_aarch64 = [-8.855500, -8.095195, -9.153704, -9.321686]
     test_list.append(cfd_flamelet_ch4_unsteady)
 
@@ -359,7 +360,8 @@ def main():
     ramp_msw.cfg_file = "inv_ramp_msw.cfg"
     ramp_msw.test_iter = 100
     ramp_msw.test_vals = [-7.242600, -1.471141, -0.077507, 0.054419]
-    ramp_msw.tol = [0.2, 0.2, 0.00001, 0.00001]
+    # Newton-Krylov residual histories vary across MPI environments; keep force checks tight.
+    ramp_msw.tol = [0.5, 0.5, 0.00001, 0.00001]
     test_list.append(ramp_msw)
 
     # MFR_coupling
@@ -485,7 +487,7 @@ def main():
     rae2822_sa.cfg_dir   = "rans/rae2822"
     rae2822_sa.cfg_file  = "turb_SA_RAE2822.cfg"
     rae2822_sa.test_iter = 20
-    rae2822_sa.test_vals = [-3.030366, -5.142544, 0.790275, 0.019666, 1000.000000]
+    rae2822_sa.test_vals = [-3.029664, -5.297480, 0.789729, 0.019682, 1000.000000]
     test_list.append(rae2822_sa)
 
     # RAE2822 SST
@@ -493,7 +495,7 @@ def main():
     rae2822_sst.cfg_dir   = "rans/rae2822"
     rae2822_sst.cfg_file  = "turb_SST_RAE2822.cfg"
     rae2822_sst.test_iter = 20
-    rae2822_sst.test_vals = [-1.485957, 5.891491, 0.675200, 0.023427, 100.000000]
+    rae2822_sst.test_vals = [-1.486489, 5.881103, 0.661483, 0.023767, 100.000000]
     test_list.append(rae2822_sst)
 
     # RAE2822 SST_SUST
@@ -593,7 +595,7 @@ def main():
     turb_naca0012_sa.cfg_dir   = "rans/naca0012"
     turb_naca0012_sa.cfg_file  = "turb_NACA0012_sa.cfg"
     turb_naca0012_sa.test_iter = 5
-    turb_naca0012_sa.test_vals = [-12.037552, -16.376949, 1.080346, 0.018385, 20.000000, -1.564082, 20.000000, -4.180917, 0.000000]
+    turb_naca0012_sa.test_vals = [-12.037577, -16.376949, 1.080346, 0.018385, 20.000000, -1.564038, 20.000000, -4.180921, 0.000000]
     turb_naca0012_sa.test_vals_aarch64 = [-12.037489, -16.376949, 1.080346, 0.018385, 20.000000, -1.564143, 20.000000, -4.180945, 0.000000]
     turb_naca0012_sa.timeout   = 3200
     test_list.append(turb_naca0012_sa)
@@ -603,17 +605,34 @@ def main():
     turb_naca0012_sst.cfg_dir   = "rans/naca0012"
     turb_naca0012_sst.cfg_file  = "turb_NACA0012_sst.cfg"
     turb_naca0012_sst.test_iter = 10
-    turb_naca0012_sst.test_vals = [-12.094692, -15.251093, -5.906365, 1.070413, 0.015775, -2.376043, 0.000000]
+    turb_naca0012_sst.test_vals = [-12.094729, -15.251094, -5.906365, 1.070413, 0.015775, -2.376188, 0.000000]
     turb_naca0012_sst.test_vals_aarch64 = [-12.075620, -15.246688, -5.861276, 1.070036, 0.015841, -1.991001, 0.000000]
     turb_naca0012_sst.timeout   = 3200
     test_list.append(turb_naca0012_sst)
 
+    # NACA0012 (SST, MUSCL_TURB with edge limiting that remains active after LIMITER_ITER).
+    turb_naca0012_sst_muscl_edge = TestCase('turb_naca0012_sst_muscl_edge')
+    turb_naca0012_sst_muscl_edge.cfg_dir = "rans/naca0012"
+    turb_naca0012_sst_muscl_edge.cfg_file = "turb_NACA0012_sst_muscl_edge.cfg"
+    turb_naca0012_sst_muscl_edge.test_iter = 10
+    turb_naca0012_sst_muscl_edge.test_vals = [-7.395492, -9.461721, -2.848756, 1.074672, 0.014912, -2.784876, 0.000000]
+    turb_naca0012_sst_muscl_edge.timeout = 3200
+    test_list.append(turb_naca0012_sst_muscl_edge)
+
+    # Keep applying the stored turbulence limiter after LIMITER_ITER.
+    turb_naca0012_sst_muscl_frozen = TestCase('turb_naca0012_sst_muscl_frozen')
+    turb_naca0012_sst_muscl_frozen.cfg_dir = "rans/naca0012"
+    turb_naca0012_sst_muscl_frozen.cfg_file = "turb_NACA0012_sst_muscl_frozen.cfg"
+    turb_naca0012_sst_muscl_frozen.test_iter = 10
+    turb_naca0012_sst_muscl_frozen.test_vals = [-8.102153, -9.310635, -3.623485, 1.071750, 0.015836, -2.241139, 0.000000]
+    test_list.append(turb_naca0012_sst_muscl_frozen)
+
     # NACA0012 (SST_SUST, FUN3D finest grid results: CL=1.0840, CD=0.01253)
     turb_naca0012_sst_sust           = TestCase('turb_naca0012_sst_sust')
     turb_naca0012_sst_sust.cfg_dir   = "rans/naca0012"
     turb_naca0012_sst_sust.cfg_file  = "turb_NACA0012_sst_sust.cfg"
     turb_naca0012_sst_sust.test_iter = 10
-    turb_naca0012_sst_sust.test_vals = [-12.081966, -14.837177, -5.733436, 1.000893, 0.019109, -2.241123]
+    turb_naca0012_sst_sust.test_vals = [-12.082100, -14.837177, -5.733436, 1.000893, 0.019109, -2.240955]
     turb_naca0012_sst_sust.test_vals_aarch64 = [-12.073964, -14.836726, -5.732390, 1.000050, 0.019144, -2.229074]
     turb_naca0012_sst_sust.timeout   = 3200
     test_list.append(turb_naca0012_sst_sust)
@@ -979,7 +998,7 @@ def main():
     turbmod_sa_neg_rae2822.cfg_dir   = "turbulence_models/sa/rae2822"
     turbmod_sa_neg_rae2822.cfg_file  = "turb_SA_NEG_RAE2822.cfg"
     turbmod_sa_neg_rae2822.test_iter = 10
-    turbmod_sa_neg_rae2822.test_vals = [1.275994, 1.062588, -2.067822, 1.300662, 0.480881, 0.000000]
+    turbmod_sa_neg_rae2822.test_vals = [1.275731, 1.063097, -2.067745, 1.301949, 0.481304, 0.000000]
     turbmod_sa_neg_rae2822.test_vals_aarch64 = [-1.345593, 1.448310, 1.208721, -0.846597, 1.248410, 0.489117, 0.000000]
     test_list.append(turbmod_sa_neg_rae2822)
 
@@ -1353,7 +1372,7 @@ def main():
     Jones_tc_restart.cfg_dir   = "turbomachinery/APU_turbocharger"
     Jones_tc_restart.cfg_file  = "Jones_restart.cfg"
     Jones_tc_restart.test_iter = 5
-    Jones_tc_restart.test_vals = [-11.905382, -12.212100, -19.212143, -13.508268, -19.078783, -13.446389, 73286.000000, 73286.000000, 0.020056, 82.286000]
+    Jones_tc_restart.test_vals = [-11.905727, -12.211201, -19.212537, -13.508616, -19.078613, -13.445675, 73286.000000, 73286.000000, 0.020056, 82.286000]
     test_list.append(Jones_tc_restart)
 
     # 2D axial stage
@@ -1461,7 +1480,7 @@ def main():
     bars_SST_2D.cfg_dir   = "sliding_interface/bars_SST_2D"
     bars_SST_2D.cfg_file  = "bars.cfg"
     bars_SST_2D.test_iter = 13
-    bars_SST_2D.test_vals = [13.000000, -0.457912, -1.541047]
+    bars_SST_2D.test_vals = [13.000000, -0.464267, -1.541047]
     bars_SST_2D.multizone = True
     test_list.append(bars_SST_2D)
 
@@ -1598,7 +1617,7 @@ def main():
     solid_periodic_pins.cfg_dir   = "solid_heat_conduction/periodic_pins"
     solid_periodic_pins.cfg_file  = "configSolid.cfg"
     solid_periodic_pins.test_iter = 750
-    solid_periodic_pins.test_vals = [-15.878960, -14.569206, 300.900000, 425.320000, 5.000000, -1.672601]
+    solid_periodic_pins.test_vals = [-15.878962, -14.569206, 300.900000, 425.320000, 5.000000, -1.672653]
     solid_periodic_pins.test_vals_aarch64 = [-15.879016, -14.569206, 300.900000, 425.320000, 5.000000, -1.672666]
     test_list.append(solid_periodic_pins)
 
@@ -1661,7 +1680,7 @@ def main():
     pywrapper_turb_naca0012_sst.cfg_dir   = "rans/naca0012"
     pywrapper_turb_naca0012_sst.cfg_file  = "turb_NACA0012_sst.cfg"
     pywrapper_turb_naca0012_sst.test_iter = 10
-    pywrapper_turb_naca0012_sst.test_vals = [-12.094692, -15.251093, -5.906365, 1.070413, 0.015775, -2.376043, 0.000000]
+    pywrapper_turb_naca0012_sst.test_vals = [-12.094729, -15.251094, -5.906365, 1.070413, 0.015775, -2.376188, 0.000000]
     pywrapper_turb_naca0012_sst.test_vals_aarch64 = [-12.075620, -15.246688, -5.861276, 1.070036, 0.015841, -1.991001, 0.000000]
     pywrapper_turb_naca0012_sst.command   = TestCase.Command("mpirun -np 2", "SU2_CFD.py", "--parallel -f")
     pywrapper_turb_naca0012_sst.timeout   = 3200
@@ -1862,9 +1881,17 @@ def main():
     species2_primitiveVenturi_mixingmodel_boundedscalar.cfg_dir   = "species_transport/venturi_primitive_3species"
     species2_primitiveVenturi_mixingmodel_boundedscalar.cfg_file  = "species2_primitiveVenturi_mixingmodel_boundedscalar.cfg"
     species2_primitiveVenturi_mixingmodel_boundedscalar.test_iter = 50
-    species2_primitiveVenturi_mixingmodel_boundedscalar.test_vals = [-5.689670, -4.511504, -4.615493, -5.795205, -0.113336, -5.704986, 5.000000, -1.433752, 5.000000, -4.921374, 5.000000, -1.771015, 0.000318, 0.000318, 0.000000, 0.000000]
+    species2_primitiveVenturi_mixingmodel_boundedscalar.test_vals = [-5.689301, -4.511142, -4.615272, -5.795460, -0.113345, -5.705122, 5.000000, -1.434208, 5.000000, -4.922106, 5.000000, -1.772145, 0.000318, 0.000318, 0.000000, 0.000000]
     test_list.append(species2_primitiveVenturi_mixingmodel_boundedscalar)
 
+    # Species edge limiting with variable density and bounded scalar transport, past LIMITER_ITER.
+    species2_primitiveVenturi_muscl_edge = TestCase('species2_primitiveVenturi_muscl_edge')
+    species2_primitiveVenturi_muscl_edge.cfg_dir = "species_transport/venturi_primitive_3species"
+    species2_primitiveVenturi_muscl_edge.cfg_file = "species2_primitiveVenturi_muscl_edge.cfg"
+    species2_primitiveVenturi_muscl_edge.test_iter = 50
+    species2_primitiveVenturi_muscl_edge.test_vals = [-5.709619, -4.522629, -4.633999, -5.815865, -0.074216, -5.653580, 5.000000, -1.427224, 5.000000, -4.986199, 5.000000, -1.791168, 0.000384, 0.000384, 0.000000, 0.000000]
+    test_list.append(species2_primitiveVenturi_muscl_edge)
+
     # 2 species (1 eq) primitive venturi mixing using mixing model including viscosity, thermal conductivity and inlet markers for SA turbulence model
     species2_primitiveVenturi_mixingmodel_viscosity           = TestCase('species2_primitiveVenturi_mixingmodel_viscosity')
     species2_primitiveVenturi_mixingmodel_viscosity.cfg_dir   = "species_transport/venturi_primitive_3species"
@@ -1894,7 +1921,7 @@ def main():
     species2_primitiveVenturi_JST.cfg_dir   = "species_transport/venturi_primitive_3species"
     species2_primitiveVenturi_JST.cfg_file  = "species2_primitiveVenturi_JST.cfg"
     species2_primitiveVenturi_JST.test_iter = 50
-    species2_primitiveVenturi_JST.test_vals = [-6.035464, -7.071918, -7.201080, -1.142940, -8.348316, 10.000000, -3.223791, 10.000000, -4.435519, 0.049048, 0.014468, 0.020068, 0.014512, 25.000000]
+    species2_primitiveVenturi_JST.test_vals = [-6.035434, -7.071855, -7.201302, -1.142934, -8.347369, 10.000000, -3.223773, 10.000000, -4.438752, 0.049046, 0.014468, 0.020067, 0.014512, 25.000000]
     test_list.append(species2_primitiveVenturi_JST)
 
     # 2 species (1 eq) primitive venturi mixing using mixing model solving enthalpy equation using preconditioning + Lax-Friedrich convective scheme
@@ -1902,7 +1929,7 @@ def main():
     species2_primitiveVenturi_Lax_Friedrich.cfg_dir   = "species_transport/venturi_primitive_3species"
     species2_primitiveVenturi_Lax_Friedrich.cfg_file  = "species2_primitiveVenturi_Lax_Friedrich.cfg"
     species2_primitiveVenturi_Lax_Friedrich.test_iter = 50
-    species2_primitiveVenturi_Lax_Friedrich.test_vals = [-6.092441, -6.981653, -6.982959, -1.195023, -8.245626, 10.000000, -3.472515, 8.000000, -5.356121, 0.048943, 0.014468, 0.020007, 0.014468, 12.500000]
+    species2_primitiveVenturi_Lax_Friedrich.test_vals = [-6.092367, -6.981726, -6.982988, -1.194909, -8.244827, 10.000000, -3.472515, 8.000000, -5.358102, 0.048940, 0.014468, 0.020005, 0.014468, 12.500000]
     test_list.append(species2_primitiveVenturi_Lax_Friedrich)
 
     # 2 species (1 eq) primitive venturi mixing
@@ -1918,7 +1945,7 @@ def main():
     species_primitiveVenturi_boundedscalar.cfg_dir     = "species_transport/venturi_primitive_3species"
     species_primitiveVenturi_boundedscalar.cfg_file    = "species2_primitiveVenturi_boundedscalar.cfg"
     species_primitiveVenturi_boundedscalar.test_iter   = 50
-    species_primitiveVenturi_boundedscalar.test_vals   = [-5.537734, -4.375388, -4.475128, -5.597904, -0.870563, -5.633513, 5.000000, -1.461524, 5.000000, -4.142508, 5.000000, -1.727874, 0.000438, 0.000438, 0.000000, 0.000000]
+    species_primitiveVenturi_boundedscalar.test_vals = [-5.537693, -4.375350, -4.475070, -5.598473, -0.870550, -5.634473, 5.000000, -1.461564, 5.000000, -4.142774, 5.000000, -1.729589, 0.000437, 0.000437, 0.000000, 0.000000]
     test_list.append(species_primitiveVenturi_boundedscalar)
 
     # 2 species (1 eq) primitive venturi mixing, unsteady and pressure-based
@@ -1935,7 +1962,7 @@ def main():
     species2_primitiveVenturi_mixingmodel_TURBULENT_MARKERS.cfg_dir   = "species_transport/venturi_primitive_3species"
     species2_primitiveVenturi_mixingmodel_TURBULENT_MARKERS.cfg_file  = "species2_primitiveVenturi_mixingmodel_TURBULENT_MARKERS.cfg"
     species2_primitiveVenturi_mixingmodel_TURBULENT_MARKERS.test_iter = 50
-    species2_primitiveVenturi_mixingmodel_TURBULENT_MARKERS.test_vals = [-4.847325, -1.922154, -1.947527, -0.716546, 1.211740, -3.930346, 21.000000, -5.152522, 9.000000, -5.325373, 4.000000, -5.865914, 2.000000, 1.000000, 0.000000, 1.000000]
+    species2_primitiveVenturi_mixingmodel_TURBULENT_MARKERS.test_vals = [-4.826194, -1.902360, -1.915879, -0.701974, 1.184582, -3.920551, 21.000000, -5.222614, 9.000000, -5.255811, 4.000000, -5.956416, 2.000000, 1.000000, 0.000000, 1.000000]
     test_list.append(species2_primitiveVenturi_mixingmodel_TURBULENT_MARKERS)
 
     # 3 species (2 eq) primitive venturi mixing with inlet files.
@@ -1960,7 +1987,7 @@ def main():
     species_passive_val.cfg_dir   = "species_transport/passive_transport_validation"
     species_passive_val.cfg_file  = "passive_transport.cfg"
     species_passive_val.test_iter = 50
-    species_passive_val.test_vals = [-16.493002, -16.246291, -16.871574, -4.257599, 10.000000, -4.526973, 8.000000, -5.193350, 0.186610, 0.000000]
+    species_passive_val.test_vals = [-16.493002, -16.246291, -16.871574, -4.253999, 10.000000, -4.526973, 8.000000, -5.184872, 0.186370, 0.000000]
     species_passive_val.test_vals_aarch64 = [-16.517744, -16.282420, -16.871663, -4.257599, 10.000000, -4.278151, 8.000000, -5.193350, 0.186610, 0.000000]
     test_list.append(species_passive_val)
 
@@ -1969,7 +1996,7 @@ def main():
     species_active_transport_temp_limits.cfg_dir   = "species_transport/passive_transport_validation"
     species_active_transport_temp_limits.cfg_file  = "active_species_transport_temp_limits.cfg"
     species_active_transport_temp_limits.test_iter = 50
-    species_active_transport_temp_limits.test_vals = [-1.785041, -2.565628, 2.460433, -3.188111, 9.000000, -5.551493, 3.000000, -5.826106, 1.456438, 0.998134, 0.001475, 0.456829]
+    species_active_transport_temp_limits.test_vals = [-1.726435, -2.525080, 2.426296, -3.132922, 9.000000, -5.514656, 3.000000, -5.742320, 1.473631, 0.997946, 0.001540, 0.474146]
     species_active_transport_temp_limits.test_vals_aarch64 = [-1.785041, -2.565628, 2.460433, -3.188111, 9.000000, -5.551493, 3.000000, -5.826106, 1.456438, 0.998134, 0.001475, 0.456829]
     test_list.append(species_active_transport_temp_limits)
 
diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py
index 250f96e66c12..6b754564c9e5 100644
--- a/TestCases/parallel_regression_AD.py
+++ b/TestCases/parallel_regression_AD.py
@@ -373,7 +373,7 @@ def main():
     discadj_flamelet_ch4_hx.cfg_file         = "lam_prem_ch4_hx_ad.cfg"
     discadj_flamelet_ch4_hx.multizone        = False
     discadj_flamelet_ch4_hx.test_iter        = 10
-    discadj_flamelet_ch4_hx.test_vals        = [-9.083502, -9.033520, -9.529644, -8.442388, -15.399568, -6.170002, -18.881156]
+    discadj_flamelet_ch4_hx.test_vals = [-9.083502, -9.033520, -9.529644, -8.442386, -15.399568, -6.169865, -18.881156]
     test_list.append(discadj_flamelet_ch4_hx)
 
     # 2D planar laminar premixed flame on isothermal burner with conjugate heat transfer (restart)
diff --git a/TestCases/rans/naca0012/turb_NACA0012_sst_muscl_edge.cfg b/TestCases/rans/naca0012/turb_NACA0012_sst_muscl_edge.cfg
new file mode 100644
index 000000000000..a157883034dd
--- /dev/null
+++ b/TestCases/rans/naca0012/turb_NACA0012_sst_muscl_edge.cfg
@@ -0,0 +1,113 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: 2D NACA 0012 Airfoil Validation Case (compressible)        %
+%                   SST with MUSCL_TURB and the VAN_ALBADA_EDGE limiter         %
+%                   http://turbmodels.larc.nasa.gov/naca0012_val_sa.html       %
+% Author: Thomas D. Economon                                                   %
+% Co-author: rois1995                                                          %
+% Institution: Stanford University                                             %
+% Date: Feb 18th, 2013                                                         %
+% File Version 8.5.0 "Harrier"                                                 %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+SOLVER= RANS
+KIND_TURB_MODEL= SST
+SST_OPTIONS= V1994m
+MATH_PROBLEM= DIRECT
+RESTART_SOL= YES
+
+% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------%
+%
+MACH_NUMBER= 0.15
+AOA= 10.0
+FREESTREAM_TEMPERATURE= 300.0
+REYNOLDS_NUMBER= 6.0E6
+REYNOLDS_LENGTH= 1.0
+
+% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
+%
+REF_ORIGIN_MOMENT_X = 0.25
+REF_ORIGIN_MOMENT_Y = 0.00
+REF_ORIGIN_MOMENT_Z = 0.00
+REF_LENGTH= 1.0
+REF_AREA= 1.0
+REF_DIMENSIONALIZATION= FREESTREAM_PRESS_EQ_ONE
+
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+MARKER_HEATFLUX= ( airfoil, 0.0 )
+MARKER_FAR= ( farfield )
+MARKER_PLOTTING= ( airfoil )
+MARKER_MONITORING= ( airfoil )
+
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
+NUM_METHOD_GRAD_RECON= LEAST_SQUARES
+CFL_NUMBER= 1000.0
+MAX_DELTA_TIME= 1E10
+CFL_ADAPT= NO
+CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 )
+ITER= 99999
+
+% ----------------------- SLOPE LIMITER DEFINITION ----------------------------%
+VENKAT_LIMITER_COEFF= 0.03
+LIMITER_ITER= 3
+
+% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
+%
+LINEAR_SOLVER= FGMRES
+LINEAR_SOLVER_PREC= ILU
+LINEAR_SOLVER_ERROR= 1E-10
+LINEAR_SOLVER_ITER= 20
+
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+CONV_NUM_METHOD_FLOW= ROE
+MUSCL_FLOW= YES
+MUSCL_KAPPA_FLOW= 0.5
+SLOPE_LIMITER_FLOW= VAN_ALBADA_EDGE
+JST_SENSOR_COEFF= ( 0.5, 0.02 )
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+
+% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------%
+%
+CONV_NUM_METHOD_TURB= SCALAR_UPWIND
+MUSCL_TURB= YES
+MUSCL_KAPPA_TURB= 0.0
+SLOPE_LIMITER_TURB= VAN_ALBADA_EDGE
+TIME_DISCRE_TURB= EULER_IMPLICIT
+CFL_REDUCTION_TURB= 1.0
+
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+% Convergence field
+CONV_FIELD= RMS_DENSITY
+CONV_RESIDUAL_MINVAL= -12
+CONV_STARTITER= 10
+CONV_CAUCHY_ELEMS= 100
+CONV_CAUCHY_EPS= 1E-6
+
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+MESH_FILENAME= n0012_225-65.su2
+MESH_FORMAT= SU2
+MESH_OUT_FILENAME= mesh_out
+SOLUTION_FILENAME= solution_flow_sst
+SOLUTION_ADJ_FILENAME= solution_adj
+TABULAR_FORMAT= CSV
+CONV_FILENAME= history
+RESTART_FILENAME= restart_flow
+RESTART_ADJ_FILENAME= restart_adj
+VOLUME_FILENAME= flow
+VOLUME_ADJ_FILENAME= adjoint
+GRAD_OBJFUNC_FILENAME= of_grad
+SURFACE_FILENAME= surface_flow
+SURFACE_ADJ_FILENAME= surface_adjoint
+OUTPUT_WRT_FREQ= 10000
+SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY, RMS_TKE, RMS_DISSIPATION, LIFT, DRAG, LINSOL_RESIDUAL, TOTAL_HEATFLUX)
+OUTPUT_FILES= (RESTART, PARAVIEW, SURFACE_PARAVIEW)
diff --git a/TestCases/rans/naca0012/turb_NACA0012_sst_muscl_frozen.cfg b/TestCases/rans/naca0012/turb_NACA0012_sst_muscl_frozen.cfg
new file mode 100644
index 000000000000..0e5fc56eea6d
--- /dev/null
+++ b/TestCases/rans/naca0012/turb_NACA0012_sst_muscl_frozen.cfg
@@ -0,0 +1,113 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: 2D NACA 0012 Airfoil Validation Case (compressible)        %
+%                   SST with MUSCL_TURB and a frozen stored limiter             %
+%                   http://turbmodels.larc.nasa.gov/naca0012_val_sa.html       %
+% Author: Thomas D. Economon                                                   %
+% Co-author: rois1995                                                          %
+% Institution: Stanford University                                             %
+% Date: Feb 18th, 2013                                                         %
+% File Version 8.5.0 "Harrier"                                                 %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+SOLVER= RANS
+KIND_TURB_MODEL= SST
+SST_OPTIONS= V1994m
+MATH_PROBLEM= DIRECT
+RESTART_SOL= YES
+
+% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------%
+%
+MACH_NUMBER= 0.15
+AOA= 10.0
+FREESTREAM_TEMPERATURE= 300.0
+REYNOLDS_NUMBER= 6.0E6
+REYNOLDS_LENGTH= 1.0
+
+% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
+%
+REF_ORIGIN_MOMENT_X = 0.25
+REF_ORIGIN_MOMENT_Y = 0.00
+REF_ORIGIN_MOMENT_Z = 0.00
+REF_LENGTH= 1.0
+REF_AREA= 1.0
+REF_DIMENSIONALIZATION= FREESTREAM_PRESS_EQ_ONE
+
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+MARKER_HEATFLUX= ( airfoil, 0.0 )
+MARKER_FAR= ( farfield )
+MARKER_PLOTTING= ( airfoil )
+MARKER_MONITORING= ( airfoil )
+
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
+NUM_METHOD_GRAD_RECON= LEAST_SQUARES
+CFL_NUMBER= 1000.0
+MAX_DELTA_TIME= 1E10
+CFL_ADAPT= NO
+CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 )
+ITER= 99999
+
+% ----------------------- SLOPE LIMITER DEFINITION ----------------------------%
+VENKAT_LIMITER_COEFF= 0.03
+LIMITER_ITER= 3
+
+% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
+%
+LINEAR_SOLVER= FGMRES
+LINEAR_SOLVER_PREC= ILU
+LINEAR_SOLVER_ERROR= 1E-10
+LINEAR_SOLVER_ITER= 20
+
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+CONV_NUM_METHOD_FLOW= ROE
+MUSCL_FLOW= YES
+MUSCL_KAPPA_FLOW= 0.5
+SLOPE_LIMITER_FLOW= NONE
+JST_SENSOR_COEFF= ( 0.5, 0.02 )
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+
+% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------%
+%
+CONV_NUM_METHOD_TURB= SCALAR_UPWIND
+MUSCL_TURB= YES
+MUSCL_KAPPA_TURB= 0.0
+SLOPE_LIMITER_TURB= VENKATAKRISHNAN_WANG
+TIME_DISCRE_TURB= EULER_IMPLICIT
+CFL_REDUCTION_TURB= 1.0
+
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+% Convergence field
+CONV_FIELD= RMS_DENSITY
+CONV_RESIDUAL_MINVAL= -12
+CONV_STARTITER= 10
+CONV_CAUCHY_ELEMS= 100
+CONV_CAUCHY_EPS= 1E-6
+
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+MESH_FILENAME= n0012_225-65.su2
+MESH_FORMAT= SU2
+MESH_OUT_FILENAME= mesh_out
+SOLUTION_FILENAME= solution_flow_sst
+SOLUTION_ADJ_FILENAME= solution_adj
+TABULAR_FORMAT= CSV
+CONV_FILENAME= history
+RESTART_FILENAME= restart_flow
+RESTART_ADJ_FILENAME= restart_adj
+VOLUME_FILENAME= flow
+VOLUME_ADJ_FILENAME= adjoint
+GRAD_OBJFUNC_FILENAME= of_grad
+SURFACE_FILENAME= surface_flow
+SURFACE_ADJ_FILENAME= surface_adjoint
+OUTPUT_WRT_FREQ= 10000
+SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY, RMS_TKE, RMS_DISSIPATION, LIFT, DRAG, LINSOL_RESIDUAL, TOTAL_HEATFLUX)
+OUTPUT_FILES= (RESTART, PARAVIEW, SURFACE_PARAVIEW)
diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py
index 1edd0378ab2a..b877b9d31f8c 100755
--- a/TestCases/serial_regression.py
+++ b/TestCases/serial_regression.py
@@ -236,7 +236,7 @@ def main():
     rae2822_sa.cfg_dir   = "rans/rae2822"
     rae2822_sa.cfg_file  = "turb_SA_RAE2822.cfg"
     rae2822_sa.test_iter = 20
-    rae2822_sa.test_vals = [-3.117674, -5.153287, 0.777070, 0.019412, 1000.000000]
+    rae2822_sa.test_vals = [-3.116257, -5.313749, 0.776303, 0.019411, 1000.000000]
     test_list.append(rae2822_sa)
 
     # RAE2822 SST
@@ -244,7 +244,7 @@ def main():
     rae2822_sst.cfg_dir   = "rans/rae2822"
     rae2822_sst.cfg_file  = "turb_SST_RAE2822.cfg"
     rae2822_sst.test_iter = 20
-    rae2822_sst.test_vals = [-1.916371, -1.461911, 5.894621, 0.611793, 0.016873, 100.000000]
+    rae2822_sst.test_vals = [-1.917411, -1.464736, 5.881403, 0.625075, 0.017288, 100.000000]
     test_list.append(rae2822_sst)
 
     # RAE2822 SST_SUST
@@ -301,7 +301,7 @@ def main():
     turb_naca0012_sa.cfg_dir   = "rans/naca0012"
     turb_naca0012_sa.cfg_file  = "turb_NACA0012_sa.cfg"
     turb_naca0012_sa.test_iter = 5
-    turb_naca0012_sa.test_vals = [-12.037319, -16.384159, 1.080346, 0.018385, 20.000000, -3.456939, 20.000000, -4.641254, 0.000000]
+    turb_naca0012_sa.test_vals = [-12.037299, -16.384159, 1.080346, 0.018385, 20.000000, -3.456989, 20.000000, -4.641253, 0.000000]
     turb_naca0012_sa.test_vals_aarch64 = [-12.037297, -16.384158, 1.080346, 0.018385, 20.000000, -3.455886, 20.000000, -4.641247, 0.000000]
     turb_naca0012_sa.timeout   = 3200
     test_list.append(turb_naca0012_sa)
@@ -311,7 +311,7 @@ def main():
     turb_naca0012_sst.cfg_dir   = "rans/naca0012"
     turb_naca0012_sst.cfg_file  = "turb_NACA0012_sst.cfg"
     turb_naca0012_sst.test_iter = 10
-    turb_naca0012_sst.test_vals = [-12.094431, -15.251082, -5.906366, 1.070413, 0.015775, -3.178469, 0.000000]
+    turb_naca0012_sst.test_vals = [-12.094460, -15.251083, -5.906366, 1.070413, 0.015775, -3.178954, 0.000000]
     turb_naca0012_sst.test_vals_aarch64 = [-12.076068, -15.246740, -5.861280, 1.070036, 0.015841, -3.297854, 0.000000]
     turb_naca0012_sst.timeout   = 3200
     test_list.append(turb_naca0012_sst)
@@ -343,7 +343,7 @@ def main():
     turb_naca0012_sst_sust_restart.cfg_dir   = "rans/naca0012"
     turb_naca0012_sst_sust_restart.cfg_file  = "turb_NACA0012_sst_sust.cfg"
     turb_naca0012_sst_sust_restart.test_iter = 10
-    turb_naca0012_sst_sust_restart.test_vals = [-12.080496, -14.837169, -5.733461, 1.000893, 0.019109, -2.634008]
+    turb_naca0012_sst_sust_restart.test_vals = [-12.080468, -14.837169, -5.733461, 1.000893, 0.019109, -2.634045]
     turb_naca0012_sst_sust_restart.test_vals_aarch64 = [-12.074189, -14.836725, -5.732398, 1.000050, 0.019144, -3.315560]
     turb_naca0012_sst_sust_restart.timeout   = 3200
     test_list.append(turb_naca0012_sst_sust_restart)
@@ -966,7 +966,7 @@ def main():
     Jones_tc_restart.cfg_dir   = "turbomachinery/APU_turbocharger"
     Jones_tc_restart.cfg_file  = "Jones_restart.cfg"
     Jones_tc_restart.test_iter = 5
-    Jones_tc_restart.test_vals = [-11.905642, -12.212181, -19.212554, -13.508053, -19.069991, -13.439494, 73286.000000, 73286.000000, 0.020056, 82.286000]
+    Jones_tc_restart.test_vals = [-11.905618, -12.211364, -19.212620, -13.507997, -19.069471, -13.439286, 73286.000000, 73286.000000, 0.020056, 82.286000]
     test_list.append(Jones_tc_restart)
 
     # 2D axial stage
@@ -1086,7 +1086,7 @@ def main():
     bars_SST_2D.cfg_dir   = "sliding_interface/bars_SST_2D"
     bars_SST_2D.cfg_file  = "bars.cfg"
     bars_SST_2D.test_iter = 13
-    bars_SST_2D.test_vals = [13.000000, -0.456143, -1.541051]
+    bars_SST_2D.test_vals = [13.000000, -0.462697, -1.541051]
     bars_SST_2D.multizone = True
     test_list.append(bars_SST_2D)
 
@@ -1164,7 +1164,7 @@ def main():
     fsi_cht.cfg_dir   = "fea_fsi/stat_fsi"
     fsi_cht.cfg_file  = "config.cfg"
     fsi_cht.test_iter = 20
-    fsi_cht.test_vals = [5.000000, -5.077002, -5.379450, -9.247804, -9.320014, -9.185034, 608.350000, -0.012973, 0.000000, 30.000000]
+    fsi_cht.test_vals = [5.000000, -5.076997, -5.379435, -9.247798, -9.319892, -9.185018, 608.350000, -0.012973, 0.000000, 30.000000]
     fsi_cht.multizone = True
     test_list.append(fsi_cht)
 
@@ -1688,7 +1688,7 @@ def main():
     pywrapper_turb_naca0012_sst.cfg_dir   = "rans/naca0012"
     pywrapper_turb_naca0012_sst.cfg_file  = "turb_NACA0012_sst.cfg"
     pywrapper_turb_naca0012_sst.test_iter = 10
-    pywrapper_turb_naca0012_sst.test_vals = [-12.094431, -15.251082, -5.906366, 1.070413, 0.015775, -3.178469, 0.000000]
+    pywrapper_turb_naca0012_sst.test_vals = [-12.094460, -15.251083, -5.906366, 1.070413, 0.015775, -3.178954, 0.000000]
     pywrapper_turb_naca0012_sst.test_vals_aarch64 = [-12.076068, -15.246740, -5.861280, 1.070036, 0.015841, -3.297854, 0.000000]
     pywrapper_turb_naca0012_sst.command   =  TestCase.Command(exec = "SU2_CFD.py", param = "-f")
     pywrapper_turb_naca0012_sst.timeout   = 3200
diff --git a/TestCases/species_transport/venturi_primitive_3species/species2_primitiveVenturi_muscl_edge.cfg b/TestCases/species_transport/venturi_primitive_3species/species2_primitiveVenturi_muscl_edge.cfg
new file mode 100644
index 000000000000..129b1997bfd5
--- /dev/null
+++ b/TestCases/species_transport/venturi_primitive_3species/species2_primitiveVenturi_muscl_edge.cfg
@@ -0,0 +1,144 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: Species mixing with 2 species, i.e. 1 transport equations  %
+%                   Including mixture dependent density, viscosity, thermal    %
+%                   conductivity, with VAN_ALBADA_EDGE species limiting        %
+% Author: Cristopher Morales Ubal                                              %
+% Co-author: rois1995                                                          %
+% Institution: Eindhoven University of Technology                              %
+% Date: 2022/06/15                                                             %
+% File Version 8.5.0 "Harrier"                                                 %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+SOLVER= INC_RANS
+KIND_TURB_MODEL= SST
+%
+% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------%
+%
+INC_DENSITY_MODEL= VARIABLE
+INC_DENSITY_INIT= 1.1766
+%
+INC_VELOCITY_INIT= ( 1.00, 0.0, 0.0 )
+%
+INC_ENERGY_EQUATION= YES
+INC_TEMPERATURE_INIT= 300.0
+%
+INC_NONDIM= DIMENSIONAL
+%
+% -------------------- FLUID PROPERTIES ------------------------------------- %
+%
+FLUID_MODEL= FLUID_MIXTURE
+%
+THERMODYNAMIC_PRESSURE= 101325.0 
+%
+MOLECULAR_WEIGHT= 28.96, 16.043
+%
+SPECIFIC_HEAT_CP = 1009.39, 1009.39
+%
+CONDUCTIVITY_MODEL= CONSTANT_CONDUCTIVITY
+THERMAL_CONDUCTIVITY_CONSTANT= 0.0357, 0.0357 
+%
+PRANDTL_LAM= 0.72, 0.72
+TURBULENT_CONDUCTIVITY_MODEL= NONE
+PRANDTL_TURB= 0.90, 0.90
+% --------------------------- VISCOSITY MODEL ---------------------------------%
+%
+VISCOSITY_MODEL= CONSTANT_VISCOSITY
+%
+MU_CONSTANT= 1.716E-5, 1.716E-5
+%
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+MARKER_HEATFLUX= ( wall, 0.0 )
+MARKER_SYM= ( axis )
+%
+SPECIFIED_INLET_PROFILE= NO
+INLET_FILENAME= inlet_venturi.dat
+INC_INLET_TYPE=  VELOCITY_INLET VELOCITY_INLET
+MARKER_INLET= ( gas_inlet,       300, 1.0, 1.0,  0.0, 0.0,\
+                air_axial_inlet, 300, 1.0, 0.0, -1.0, 0.0 )
+MARKER_INLET_SPECIES= (gas_inlet, 1.0,\          
+                       air_axial_inlet, 0.6)           
+%
+INC_OUTLET_TYPE= PRESSURE_OUTLET
+MARKER_OUTLET= ( outlet, 0.0)
+%
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
+%
+CFL_NUMBER= 60
+CFL_REDUCTION_SPECIES= 1.0
+CFL_REDUCTION_TURB= 1.0
+%
+ITER= 1000
+%
+% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
+%
+LINEAR_SOLVER= FGMRES
+LINEAR_SOLVER_PREC= ILU
+LINEAR_SOLVER_ERROR= 1E-8
+LINEAR_SOLVER_ITER= 5
+%
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+CONV_NUM_METHOD_FLOW= FDS
+MUSCL_FLOW= YES
+SLOPE_LIMITER_FLOW = NONE
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+%
+% -------------------- SCALAR TRANSPORT ---------------------------------------%
+%
+KIND_SCALAR_MODEL= SPECIES_TRANSPORT
+DIFFUSIVITY_MODEL= CONSTANT_DIFFUSIVITY
+DIFFUSIVITY_CONSTANT= 0.001
+%
+CONV_NUM_METHOD_SPECIES= BOUNDED_SCALAR
+MUSCL_SPECIES= YES
+SLOPE_LIMITER_SPECIES= VAN_ALBADA_EDGE
+LIMITER_ITER= 3
+%
+TIME_DISCRE_SPECIES= EULER_IMPLICIT
+%
+SPECIES_INIT= 1.0
+SPECIES_CLIPPING= YES
+SPECIES_CLIPPING_MIN= 0.0
+SPECIES_CLIPPING_MAX= 1.0
+%
+% -------------------- TURBULENT TRANSPORT ---------------------------------------%
+%
+CONV_NUM_METHOD_TURB= BOUNDED_SCALAR
+MUSCL_TURB= NO
+%
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+CONV_FIELD= RMS_PRESSURE, RMS_VELOCITY-X, RMS_VELOCITY-Y, RMS_TKE, RMS_SPECIES
+CONV_RESIDUAL_MINVAL= -18
+CONV_STARTITER= 10
+%
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+MESH_FILENAME= primitiveVenturi.su2
+SCREEN_OUTPUT= INNER_ITER WALL_TIME \
+               RMS_PRESSURE RMS_VELOCITY-X RMS_VELOCITY-Y RMS_TKE RMS_DISSIPATION RMS_SPECIES_0 \
+               LINSOL_ITER LINSOL_RESIDUAL \
+               LINSOL_ITER_TURB LINSOL_RESIDUAL_TURB \
+               LINSOL_ITER_SPECIES LINSOL_RESIDUAL_SPECIES SURFACE_SPECIES_VARIANCE
+SCREEN_WRT_FREQ_INNER= 10
+%
+HISTORY_OUTPUT= RMS_RES FLOW_COEFF LINSOL SPECIES_COEFF SPECIES_COEFF_SURF
+MARKER_ANALYZE= outlet gas_inlet air_axial_inlet
+%
+OUTPUT_FILES= RESTART_ASCII, PARAVIEW_MULTIBLOCK
+VOLUME_OUTPUT= RESIDUAL, PRIMITIVE
+OUTPUT_WRT_FREQ= 1000
+%
+RESTART_SOL= NO
+SOLUTION_FILENAME= solution
+RESTART_FILENAME= restart
+%
+WRT_PERFORMANCE= YES
diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py
index 22e9e55fdb03..28cf7ba38cbd 100644
--- a/TestCases/tutorials.py
+++ b/TestCases/tutorials.py
@@ -165,7 +165,7 @@ def main():
     kenics_mixer_tutorial.cfg_dir   = "../Tutorials/incompressible_flow/Inc_Species_Transport_Composition_Dependent_Model"
     kenics_mixer_tutorial.cfg_file  = "kenics_mixer_tutorial.cfg"
     kenics_mixer_tutorial.test_iter = 10
-    kenics_mixer_tutorial.test_vals = [-7.490448, -6.823978, -6.838620, -6.383873, -7.878123, -3.003654, -7.452198, 5.000000, -1.857226, 4.000000, -5.336949, 3.000000, -6.363373, 0.025668, 0.000000, 0.025668, 0.000000, 62.717000, 8.462600, 46.725000, 7.529400]
+    kenics_mixer_tutorial.test_vals = [-7.490463, -6.823949, -6.838579, -6.383927, -7.879236, -3.003913, -7.452024, 5.000000, -1.857322, 4.000000, -5.342039, 3.000000, -6.366408, 0.025672, 0.000000, 0.025672, 0.000000, 62.693000, 8.462600, 46.701000, 7.529400]
     kenics_mixer_tutorial.command   = TestCase.Command("mpirun -n 2", "SU2_CFD")
     test_list.append(kenics_mixer_tutorial)
 
@@ -259,7 +259,7 @@ def main():
     tutorial_trans_flatplate_T3Am.cfg_dir    = "../Tutorials/compressible_flow/Transitional_Flat_Plate/Langtry_and_Menter/T3A-"
     tutorial_trans_flatplate_T3Am.cfg_file   = "transitional_LM_model_ConfigFile.cfg"
     tutorial_trans_flatplate_T3Am.test_iter  = 20
-    tutorial_trans_flatplate_T3Am.test_vals  = [-5.587389, -1.700868, -3.093935, -0.102834, -3.750523, 3.287643, -2.394575, 1.119623]
+    tutorial_trans_flatplate_T3Am.test_vals = [-5.587379, -1.700868, -3.093915, -0.102826, -3.750523, 3.287643, -2.394575, 1.119623]
     tutorial_trans_flatplate_T3Am.test_vals_aarch64 = [-5.540938, -1.681627, -2.878831, -0.058224, -3.695533, 3.413628, -2.385345, 1.103633]
     tutorial_trans_flatplate_T3Am.no_restart = True
     test_list.append(tutorial_trans_flatplate_T3Am)
@@ -314,7 +314,7 @@ def main():
     tutorial_unst_naca0012.cfg_dir       = "../Tutorials/compressible_flow/Unsteady_NACA0012"
     tutorial_unst_naca0012.cfg_file      = "unsteady_naca0012.cfg"
     tutorial_unst_naca0012.test_iter     = 520
-    tutorial_unst_naca0012.test_vals = [520.000000, 0.000000, -5.292186, 0.000000, 0.304255, 0.792893, 0.001907, 0.010888]
+    tutorial_unst_naca0012.test_vals = [520.000000, 0.000000, -5.292091, 0.000000, 0.312095, 0.781939, 0.001784, 0.010509]
     tutorial_unst_naca0012.test_vals_aarch64 = [520, 0, -5.292359, 0, 0.284720, 0.766329, 0.000954, 0.007565]
     tutorial_unst_naca0012.unsteady      = True
     test_list.append(tutorial_unst_naca0012)
diff --git a/TestCases/vandv.py b/TestCases/vandv.py
index 88602b2d2275..b5ade77f7497 100644
--- a/TestCases/vandv.py
+++ b/TestCases/vandv.py
@@ -63,7 +63,7 @@ def main():
     flatplate_sst1994m.cfg_dir   = "vandv/rans/flatplate"
     flatplate_sst1994m.cfg_file  = "turb_flatplate_sst.cfg"
     flatplate_sst1994m.test_iter = 5
-    flatplate_sst1994m.test_vals = [-13.040943, -10.136971, -10.942764, -7.985181, -10.323857, -4.732487, 0.002801]
+    flatplate_sst1994m.test_vals = [-13.041030, -10.136996, -10.941711, -7.986212, -10.323862, -4.732531, 0.002801]
     flatplate_sst1994m.test_vals_aarch64 = [-13.021715, -9.534786, -10.401912, -7.501836, -9.750800, -4.850665, 0.002807]
     test_list.append(flatplate_sst1994m)
 
@@ -72,7 +72,7 @@ def main():
     bump_sst1994m.cfg_dir   = "vandv/rans/bump_in_channel"
     bump_sst1994m.cfg_file  = "turb_bump_sst.cfg"
     bump_sst1994m.test_iter = 5
-    bump_sst1994m.test_vals = [-11.928292, -10.095796, -9.512953, -6.445652, -11.774088, -6.988752, 0.004931]
+    bump_sst1994m.test_vals = [-11.928276, -10.095787, -9.512963, -6.445637, -11.774253, -6.985525, 0.004931]
     bump_sst1994m.test_vals_aarch64 = [-13.042689, -10.812982, -10.604523, -7.655547, -10.816257, -5.308083, 0.004911]
     test_list.append(bump_sst1994m)
 
@@ -99,7 +99,7 @@ def main():
     dsma661_sa.cfg_dir    = "vandv/rans/dsma661"
     dsma661_sa.cfg_file   = "dsma661_sa_config.cfg"
     dsma661_sa.test_iter  = 5
-    dsma661_sa.test_vals  = [-11.255379, -8.243275, -9.037314, -5.927070, -10.737675, 0.155687, 0.024232]
+    dsma661_sa.test_vals = [-11.255401, -8.243277, -9.037333, -5.927120, -10.737675, 0.155687, 0.024232]
     dsma661_sa.test_vals_aarch64 = [-11.293183, -8.241775, -9.083761, -6.011398, -10.737680, 0.155687, 0.024232]
     test_list.append(dsma661_sa)
 
@@ -108,7 +108,7 @@ def main():
     dsma661_sst.cfg_dir   = "vandv/rans/dsma661"
     dsma661_sst.cfg_file  = "dsma661_sst_config.cfg"
     dsma661_sst.test_iter = 5
-    dsma661_sst.test_vals = [-11.025153, -8.156995, -9.057021, -5.947228, -10.650874, -7.884423, 0.155882, 0.023344]
+    dsma661_sst.test_vals = [-11.025513, -8.156985, -9.056477, -5.949486, -10.650930, -7.888614, 0.155882, 0.023344]
     dsma661_sst.test_vals_aarch64 = [-10.977195, -8.403731, -8.747068, -5.808899, -10.522786, -7.369851, 0.155875, 0.023353]
     test_list.append(dsma661_sst)
 
@@ -121,7 +121,7 @@ def main():
     sandiajet_sst.cfg_dir   = "vandv/species_transport/sandia_jet"
     sandiajet_sst.cfg_file  = "validation.cfg"
     sandiajet_sst.test_iter = 5
-    sandiajet_sst.test_vals = [-4.930056, -3.072757, -3.531917, -3.632319, 1.186060, -6.175518, 5.000000, -1.413277, 5.000000, -4.735257, 5.000000, -1.924560, 0.000256, 0.000000, 0.000000, 0.000256, 4126.200000, 4008.100000, 62.399000, 55.713000]
+    sandiajet_sst.test_vals = [-4.930056, -3.072757, -3.531918, -3.632508, 1.186029, -6.175346, 5.000000, -1.413277, 5.000000, -4.734818, 5.000000, -1.923974, 0.000256, 0.000000, 0.000000, 0.000256, 4126.200000, 4008.100000, 62.399000, 55.713000]
     sandiajet_sst.test_vals_aarch64 = [-8.271400, -6.305945, -6.123143, -5.634565, -1.687556, -7.869775, 5.000000, -1.527449, 5.000000, -4.053788, 5.000000, -2.303502, 0.000257, 0.000000, 0.000000, 0.000257, 4020.500000, 3919.900000, 49.151000, 51.436000]
     test_list.append(sandiajet_sst)
 
diff --git a/UnitTests/SU2_CFD/scalar_muscl.cpp b/UnitTests/SU2_CFD/scalar_muscl.cpp
new file mode 100644
index 000000000000..7180bbdb3be3
--- /dev/null
+++ b/UnitTests/SU2_CFD/scalar_muscl.cpp
@@ -0,0 +1,163 @@
+/*!
+ * \file scalar_muscl.cpp
+ * \brief Tests for scalar reconstruction, boundary fluxes and gradient storage.
+ * \author rois1995
+ * \version 8.5.0 "Harrier"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md)
+ *
+ * SU2 is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * SU2 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "catch.hpp"
+#include "../UnitQuadTestCase.hpp"
+#include "../../SU2_CFD/include/solvers/CTurbSolver.hpp"
+#include "../../SU2_CFD/include/numerics/scalar/scalar_edge_flux.hpp"
+#include "../../SU2_CFD/include/variables/CSpeciesVariable.hpp"
+#include "../../SU2_CFD/include/variables/CHeatVariable.hpp"
+#include "../../SU2_CFD/include/variables/CPoissonVariable.hpp"
+
+TEST_CASE("Edge reconstruction of scalars and of the flow", "[ScalarMuscl]") {
+  auto edgeFace = [](auto value, auto edgeLimiter, su2double delta, su2double projection) {
+    using Double = decltype(value);
+    using Int = typename CLaneTraits<Double>::Int;
+    C3DDoubleMatrix gradient(2, 1, 2);
+    gradient(0, 0, 0) = projection;
+    gradient(1, 0, 0) = -projection;
+    Vector<Double, 2> edge;
+    edge(0) = Double(1.0);
+    edge(1) = Double(0.0);
+    CPair<CScalarValues<Double, 1>> face;
+    face.i.all(0) = Double(0.0);
+    face.j.all(0) = Double(delta);
+    su2activematrix unusedLimiter(2, 1);
+    reconstruct<1, decltype(edgeLimiter)::value>(Int(0UL), Int(1UL), edge, gradient, unusedLimiter,
+                                                 LIMITER::VAN_ALBADA_EDGE, 0, face, Double(0.0), Double(1.0));
+    return face;
+  };
+  /*--- Scalars use the van Albada function of the point limiters, the flow keeps its own form. ---*/
+  const std::integral_constant<MusclLimiter, MusclLimiter::EDGE_SCALAR> scalarEdge{};
+  const std::integral_constant<MusclLimiter, MusclLimiter::EDGE> flowEdge{};
+  auto scalarLimiter = [](su2double proj, su2double delta) {
+    return LimiterHelpers<>::vanAlbadaFunction(proj, delta, EPS);
+  };
+  auto flowLimiter = [](su2double proj, su2double delta) {
+    const su2double delta2 = pow(delta, 2) + 1e-6;
+    return (delta2 + proj * delta) / (pow(proj, 2) + delta2);
+  };
+  /*--- Some faces cancel to (nearly) zero, where only rounding is left (it depends on the compiler, e.g. on
+   *    FMA contraction). The margin is far above that rounding and far below the other face values (>= 1e-9). ---*/
+  const auto approx = [](const auto& value) { return Approx(value).epsilon(1e-12).margin(1e-15); };
+  for (const su2double delta : {0.0, 1e-8, -1e-8, 1e-3, -1e-3}) {
+    for (const su2double projection : {0.0, 2e-8, -2e-8, 2e-3}) {
+      const auto checkFace = [&](auto edgeLimiter, auto limiter) {
+        const auto scalar = edgeFace(su2double(0.0), edgeLimiter, delta, projection);
+        const auto vector = edgeFace(simd::Array<su2double, 4>(0.0), edgeLimiter, delta, projection);
+        const auto left = 0.5 * projection * limiter(projection, delta);
+        const auto right = delta + 0.5 * projection * limiter(-projection, delta);
+        CHECK(scalar.i.all(0) == approx(left));
+        CHECK(scalar.j.all(0) == approx(right));
+        for (size_t lane = 0; lane < 4; ++lane) {
+          CHECK(vector.i.all(0)[lane] == approx(left));
+          CHECK(vector.j.all(0)[lane] == approx(right));
+        }
+      };
+      checkFace(scalarEdge, scalarLimiter);
+      checkFace(flowEdge, flowLimiter);
+    }
+  }
+}
+
+struct ScalarBoundaryTestSolver : CScalarSolver<CTurbVariable> {
+  using CScalarSolver::BoundedScalarBCFlux;
+  void LoadRestart(CGeometry**, CSolver***, CConfig*, int, bool) override {}
+  ScalarBoundaryTestSolver(CGeometry* geometry, CConfig* config, const CSolver* flow, bool conservative)
+      : CScalarSolver(geometry, config, flow, conservative, true) {
+    nPoint = geometry->GetnPoint();
+    nPointDomain = geometry->GetnPointDomain();
+    nVar = 1;
+    nodes = new CTurbVariable(nPoint, nDim, nVar, config);
+    SetBaseClassPointerToNodes();
+    LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0);
+    Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config);
+  }
+};
+
+TEST_CASE("Bounded scalar boundary flux upwinds density and weights the correction", "[ScalarMuscl]") {
+  UnitQuadTestCase test;
+  test.InitConfig();
+  test.InitGeometry();
+  test.InitSolver();
+  for (const bool conservative : {false, true}) {
+    ScalarBoundaryTestSolver solver(test.geometry.get(), test.config.get(), test.solver[FLOW_SOL], conservative);
+    solver.GetNodes()->GetSolution(0)[0] = 0.5;
+    const su2double rho = 2.0, ghostRho = 5.0, normal[3] = {1.0, 0.0, 0.0};
+    for (const bool implicit : {false, true}) {
+      for (const su2double speed : {-3.0, 0.0, 3.0}) {
+        for (const su2double weight : {0.25, 1.0}) {
+          for (const bool ghost : {false, true}) {
+            solver.LinSysRes.SetValZero();
+            solver.Jacobian.SetValZero();
+            const su2double velocity[3] = {speed, 0.0, 0.0};
+            const auto mass =
+                solver.BoundedScalarBCFlux(0, implicit, rho, velocity, normal, ghost ? &ghostRho : nullptr, weight);
+            const auto expectedMass = speed * (ghost && speed < 0.0 ? ghostRho : rho);
+            const auto q = weight * (conservative ? expectedMass : speed);
+            CHECK(mass == Approx(expectedMass));
+            CHECK(solver.LinSysRes(0, 0) == Approx(-0.5 * q));
+            CHECK(solver.Jacobian.GetBlockView(0, 0)(0, 0) == Approx(implicit ? -q / (conservative ? rho : 1.0) : 0.0));
+          }
+        }
+      }
+    }
+  }
+}
+
+TEST_CASE("Scalar reconstruction gradients follow the constructor's solver-specific flag", "[ScalarMuscl]") {
+  struct Variable : CScalarVariable {
+    using CScalarVariable::CScalarVariable;
+    bool HasAuxiliaryGradient() const { return Gradient_Aux.size() != 0; }
+  };
+  for (const bool separate : {false, true}) {
+    for (const bool scalarMuscl : {false, true}) {
+      std::stringstream options;
+      options << "SOLVER= EULER\nMUSCL_FLOW= YES\nNUM_METHOD_GRAD= GREEN_GAUSS\n"
+              << "MUSCL_TURB= " << (scalarMuscl ? "YES" : "NO") << "\n"
+              << "MUSCL_SPECIES= " << (scalarMuscl ? "YES" : "NO") << "\n"
+              << "MUSCL_HEAT= " << (scalarMuscl ? "YES" : "NO") << "\n"
+              << "NUM_METHOD_GRAD_RECON= " << (separate ? "LEAST_SQUARES" : "GREEN_GAUSS") << "\n";
+      CConfig config(options, SU2_COMPONENT::SU2_CFD, false);
+      REQUIRE(config.GetReconstructionGradientRequired() == separate);
+      for (const bool muscl : {false, true}) {
+        Variable variable(2, 2, 1, &config, muscl);
+        CHECK(variable.HasAuxiliaryGradient() == (separate && muscl));
+        CHECK((&variable.GetGradient_Reconstruction() != &variable.GetGradient()) == (separate && muscl));
+      }
+      CTurbVariable turb(2, 2, 1, &config);
+      const su2double speciesInf[1] = {0.5};
+      CSpeciesVariable species(speciesInf, 2, 2, 1, &config);
+      CHeatVariable heat(1.0, 2, 2, 1, &config);
+      CPoissonVariable poisson(0.0, 2, 2, 1, &config);
+      CHECK((&turb.GetGradient_Reconstruction() != &turb.GetGradient()) == (separate && scalarMuscl));
+      CHECK((&species.GetGradient_Reconstruction() != &species.GetGradient()) == (separate && scalarMuscl));
+      CHECK((&heat.GetGradient_Reconstruction() != &heat.GetGradient()) == (separate && scalarMuscl));
+      CHECK(&poisson.GetGradient_Reconstruction() == &poisson.GetGradient());
+    }
+  }
+}
diff --git a/UnitTests/UnitQuadTestCase.hpp b/UnitTests/UnitQuadTestCase.hpp
index d8da75dc0960..cb8a4aa56b68 100644
--- a/UnitTests/UnitQuadTestCase.hpp
+++ b/UnitTests/UnitQuadTestCase.hpp
@@ -87,6 +87,7 @@ struct UnitQuadTestCase {
     cout.rdbuf(nullptr);
     {
       auto aux_geometry = std::unique_ptr<CGeometry>(new CPhysicalGeometry(config.get(), 0, 1));
+      aux_geometry->SetColorGrid_Parallel(config.get());
       geometry = std::unique_ptr<CGeometry>(new CPhysicalGeometry(aux_geometry.get(), config.get()));
     }
     geometry->SetSendReceive(config.get());
diff --git a/UnitTests/meson.build b/UnitTests/meson.build
index 2c9acd98b03d..1883237365ec 100644
--- a/UnitTests/meson.build
+++ b/UnitTests/meson.build
@@ -21,6 +21,7 @@ su2_cfd_tests = files(['Common/CConfig_tests.cpp',
                        'SU2_CFD/output/COutput_convergence_tests.cpp',
                        'SU2_CFD/gradients.cpp',
                        'SU2_CFD/periodic_limiters.cpp',
+                       'SU2_CFD/scalar_muscl.cpp',
                        'SU2_CFD/solvers/CNEMOEulerSolver_tests.cpp',
                        'SU2_CFD/nemo_viscous_assembly.cpp',
                        'SU2_CFD/windowing.cpp',
diff --git a/config_template.cfg b/config_template.cfg
index 0aa29398b90a..b44fd8797a0d 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -1663,7 +1663,7 @@ MUSCL_TURB= NO
 %
 MUSCL_KAPPA_TURB= 0.0
 %
-% Slope limiter (same as SLOPE_LIMITER_FLOW except VAN_ALBADA_EDGE)
+% Slope limiter (same as SLOPE_LIMITER_FLOW)
 %
 SLOPE_LIMITER_TURB= VENKATAKRISHNAN
 %
@@ -1703,7 +1703,7 @@ SENS_REMOVE_SHARP = NO
 %
 SENS_SMOOTHING= NONE
 %
-% Freeze the value of the limiter after a number of iterations
+% Freeze stored point limiters after this iteration; VAN_ALBADA_EDGE is always recomputed.
 LIMITER_ITER= 999999
 %
 % 1st order artificial dissipation coefficients for