diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index dc96c21ad3b7..d2435111601d 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -1097,6 +1097,19 @@ class CConfig {
   su2double Const_DES;                 /*!< \brief Detached Eddy Simulation Constant. */
   WINDOW_FUNCTION Kind_WindowFct;      /*!< \brief Type of window (weight) function for objective functional. */
   unsigned short Kind_HybridRANSLES;   /*!< \brief Kind of Hybrid RANS/LES. */
+  bool StochasticBackscatter;          /*!< \brief Option to include Stochastic Backscatter Model. */
+  su2double SBS_Cdelta;                /*!< \brief Stochastic Backscatter Model lengthscale coefficient. */
+  unsigned short SBS_maxIterSmooth;    /*!< \brief Maximum number of smoothing iterations for the SBS model. */
+  su2double SBS_Ctau;                  /*!< \brief Stochastic Backscatter Model timescale coefficient. */
+  su2double SBS_Cmag;                  /*!< \brief Stochastic Backscatter Model intensity coefficient. */
+  bool stochSourceNu;                  /*!< \brief Option for including stochastic source term in turbulence model equation (Stochastic Backscatter Model). */
+  bool stochSourceDiagnostics;         /*!< \brief Option for writing diagnostics related to stochastic source terms in Langevin equations (Stochastic Backscatter Model). */
+  bool StochBackscatterInBox;          /*!< \brief Option for activating the Stochastic Backscatter Model only in a bounded box. */
+  su2double StochBackscatterBoxBounds[6]; /*!< \brief Bounds of the box where the Stochastic Backscatter Model is active. */
+  su2double stochFdThreshold;          /*!< \brief Shielding function lower threshold for application of Stochastic Backscatter Model. */
+  su2double stochSourceRelax;          /*!< \brief Relaxation factor for stochastic source term generation (Stochastic Backscatter Model). */
+  bool enforceLES;                     /*!< \brief Option to enforce LES mode in hybrid RANS-LES simulations. */
+  su2double LES_FilterWidth;           /*!< \brief LES filter width for hybrid RANS-LES simulations. */
   unsigned short Kind_RoeLowDiss;      /*!< \brief Kind of Roe scheme with low dissipation for unsteady flows. */
 
   unsigned short nSpanWiseSections; /*!< \brief number of span-wise sections */
@@ -9595,6 +9608,42 @@ class CConfig {
    */
   unsigned short GetKind_HybridRANSLES(void) const { return Kind_HybridRANSLES; }
 
+  /*!
+   * \brief Get if the Stochastic Backscatter Model must be activated.
+   * \return TRUE if the Stochastic Backscatter Model is activated.
+   */
+  bool GetStochastic_Backscatter(void) const { return StochasticBackscatter; }
+
+  /*!
+   * \brief Get if the LES mode must be enforced.
+   * \return TRUE if LES is enforced.
+   */
+  bool GetEnforceLES(void) const { return enforceLES; }
+
+  /*!
+   * \brief Get if the stochastic source term must be included in the turbulence model equation.
+   * \return TRUE if the stochastic source term is included in the turbulence model equation.
+   */
+  bool GetStochSourceNu(void) const { return stochSourceNu; }
+
+  /*!
+   * \brief Get if the diagnostics of the stochastic source terms in Langevin equations must be computed.
+   * \return TRUE if the diagnostics of the stochastic source terms in Langevin equations are computed.
+   */
+  bool GetStochSourceDiagnostics(void) const { return stochSourceDiagnostics; }
+
+  /*!
+   * \brief Get the relaxation factor for stochastic source term generation.
+   * \return Relaxation factor for stochastic source term generation.
+   */
+  su2double GetStochSourceRelax(void) const { return stochSourceRelax; }
+
+  /*!
+   * \brief Get the LES Filter Width.
+   * \return Value of LES Filter Width.
+   */
+  su2double GetLES_FilterWidth(void) const { return LES_FilterWidth; }
+
   /*!
    * \brief Get the Kind of Roe Low Dissipation Scheme for Unsteady flows.
    * \return Value of Low dissipation approach.
@@ -9607,6 +9656,48 @@ class CConfig {
    */
   su2double GetConst_DES(void) const { return Const_DES; }
 
+  /*!
+   * \brief Get the SBS lengthscale coefficient.
+   * \return Value of SBS lengthscale coefficient.
+   */
+  su2double GetSBS_Cdelta(void) const { return SBS_Cdelta; }
+
+  /*!
+   * \brief Get the SBS timescale coefficient.
+   * \return Value of SBS timescale coefficient.
+   */
+  unsigned short GetSBS_maxIterSmooth(void) const { return SBS_maxIterSmooth; }
+
+  /*!
+   * \brief Get the SBS timescale coefficient.
+   * \return Value of SBS timescale coefficient.
+   */
+  su2double GetSBS_Ctau(void) const { return SBS_Ctau; }
+
+  /*!
+   * \brief Get the SBS intensity coefficient.
+   * \return Value of SBS intensity coefficient.
+   */
+  su2double GetSBS_Cmag(void) const { return SBS_Cmag; }
+
+  /*!
+   * \brief Get if the Stochastic Backscatter Model must be applied only in a bounded box.
+   * \return True if the model is applied in a bounded box.
+   */
+  bool GetStochBackscatterInBox(void) const { return StochBackscatterInBox; }
+
+  /*!
+   * \brief Get the bounds of the box where the Stochastic Backscatter Model is active.
+   * \return Bounds of the box where the model is active.
+   */
+  const su2double* GetStochBackscatterBoxBounds(void) const { return StochBackscatterBoxBounds; }
+
+  /*!
+   * \brief Get shielding function lower threshold for application of Stochastic Backscatter Model.
+   * \return Shielding function threshold for application of the model.
+   */
+  su2double GetStochFdThreshold(void) const { return stochFdThreshold; }
+
   /*!
    * \brief Get the type of tape that will be checked in a tape debug run.
    */
diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index 690350fae050..548871ba9d17 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -820,7 +820,8 @@ enum class CENTERED {
   JST,            /*!< \brief Jameson-Smith-Turkel centered numerical method. */
   LAX,            /*!< \brief Lax-Friedrich centered numerical method. */
   JST_MAT,        /*!< \brief JST with matrix dissipation. */
-  JST_KE          /*!< \brief Kinetic Energy preserving Jameson-Smith-Turkel centered numerical method. */
+  JST_KE,         /*!< \brief Kinetic Energy preserving Jameson-Smith-Turkel centered numerical method. */
+  LD2             /*!< \brief Low-Dissipation Low-Dispersion (LD2) centered scheme. */
 };
 static const MapType<std::string, CENTERED> Centered_Map = {
   MakePair("NONE", CENTERED::NONE)
@@ -828,6 +829,7 @@ static const MapType<std::string, CENTERED> Centered_Map = {
   MakePair("JST_KE", CENTERED::JST_KE)
   MakePair("JST_MAT", CENTERED::JST_MAT)
   MakePair("LAX-FRIEDRICH", CENTERED::LAX)
+  MakePair("LD2", CENTERED::LD2)
 };
 
 
@@ -2697,6 +2699,7 @@ enum class MPI_QUANTITIES {
   MAX_LENGTH           ,  /*!< \brief Maximum length communication. */
   GRID_VELOCITY        ,  /*!< \brief Grid velocity communication. */
   SOLUTION_EDDY        ,  /*!< \brief Turbulent solution plus eddy viscosity communication. */
+  STOCH_SOURCE_LANG    ,  /*!< \brief Stochastic source term for Langevin equations communication. */
   SOLUTION_MATRIX      ,  /*!< \brief Matrix solution communication. */
   SOLUTION_MATRIXTRANS ,  /*!< \brief Matrix transposed solution communication. */
   NEIGHBORS            ,  /*!< \brief Neighbor point count communication (for JST). */
diff --git a/Common/include/toolboxes/random_toolbox.hpp b/Common/include/toolboxes/random_toolbox.hpp
new file mode 100644
index 000000000000..0249e1c5076f
--- /dev/null
+++ b/Common/include/toolboxes/random_toolbox.hpp
@@ -0,0 +1,153 @@
+/*!
+ * \file random_toolbox.hpp
+ * \brief Collection of utility functions for random number generation.
+ * \version 8.3.0 "Harrier"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2025, 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/>.
+ */
+
+#pragma once
+#include <cstdint>
+
+namespace RandomToolbox {
+/// \addtogroup RandomToolbox
+/// @{
+
+/*!
+ * \brief SplitMix64 hash function for 64-bit integers.
+ * \param[in] x Input value to hash.
+ * \return Hashed 64-bit output.
+ */
+static inline uint64_t splitmix64(uint64_t x) {
+    x += 0x9e3779b97f4a7c15ULL;                  // golden ratio offset
+    x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL; // first mixing step
+    x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL; // second mixing step
+    return x ^ (x >> 31);                        // final avalanche
+}
+
+/*!
+ * \brief Generate a deterministic 64-bit hash from three integers.
+ * \param[in] nodeIndex Global node index.
+ * \param[in] iDim Dimension index.
+ * \param[in] timeIter Current time iteration of the simulation.
+ * \return 64-bit hash value.
+ */
+inline uint64_t GetHash(unsigned long nodeIndex, unsigned short iDim, unsigned long timeIter) {
+    uint64_t x = nodeIndex;
+    x ^= splitmix64(iDim);
+    x ^= splitmix64(timeIter);
+    return splitmix64(x);
+}
+
+/*!
+ * \brief Convert a 64-bit hash into a uniform double in (0,1].
+ * Uses the top 53 bits of the hash to fill the mantissa of a double.
+ * Ensures the result is never zero, suitable for Box-Muller transform.
+ * \param[in] x 64-bit hash.
+ * \return Uniform double in the interval (0,1].
+ */
+inline double HashToUniform(uint64_t x) {
+    constexpr double inv53 = 1.0 / 9007199254740992.0; // 1/2^53
+    uint64_t uInt = x >> 11;                           // top 53 bits
+    return (uInt + 1) * inv53;                         // map to (0,1]
+}
+
+/*!
+ * \brief Generate a standard normal random number from a 64-bit hash.
+ * Uses two deterministic uniforms derived from the hash and its bitwise NOT
+ * as inputs to the Box-Muller transform.
+ * \param[in] x 64-bit hash.
+ * \return Standard normal random number (mean=0, stddev=1).
+ */
+inline double HashToNormal(uint64_t x) {
+    constexpr double pi = 3.14159265358979323846;
+    double u = HashToUniform(x);    // first uniform
+    double v = HashToUniform(~x);   // second uniform (bitwise NOT)
+    double r = sqrt(-2.0 * log(u));  
+    double theta = 2.0 * pi * v;
+    return r * cos(theta);          // one normal sample
+}
+
+/*!
+ * \brief Generate a deterministic standard normal number for a cell, dimension, and timestep.
+ * 
+ * Combines hashing and Box-Muller in one function.
+ * 
+ * \param[in] nodeIndex Global node index.
+ * \param[in] dim Dimension index.
+ * \param[in] timeIter Simulation timestep (1-based).
+ * \return Standard normal random number.
+ */
+inline double GetNormal(unsigned long nodeIndex, unsigned long dim, unsigned long timeIter) {
+    uint64_t hash = GetHash(nodeIndex, dim, timeIter);
+    return HashToNormal(hash);
+}
+
+/*!
+ * \brief Compute modified bessel function of first kind (order 0).
+ * \param[in] x Argument of Bessel funtion.
+ * \return Value of Bessel function.
+ */
+inline double GetBesselZero(double x) {
+  double abx = fabs(x);
+  if (abx < 3.75) {
+    double t = abx / 3.75;
+    double p = 1.0 + t * t * (3.5156229 + t * t * (3.0899424 + t * t * (1.2067492 + t * t * (0.2659732 + t * t * (0.0360768 + t * t * 0.0045813)))));
+    return log(p);
+  } else {
+    double t = 3.75 / abx;
+    double poly = 0.39894228 + t * (0.01328592 + t * (0.00225319 + t * (-0.00157565 + t * (0.00916281 + t * (-0.02057706 + t * (0.02635537 + t * (-0.01647633 + t * 0.00392377)))))));
+    return abx - 0.5 * log(abx) + log(poly);
+  }
+}
+
+/*!
+ * \brief Compute integral involving the product of three modified Bessel functions. 
+ *  Useful for scaling the smoothed stochastic source terms in Langevin equations.
+ * \param[in] beta_x Argument in x-direction.
+ * \param[in] beta_y Argument in y-direction.
+ * \param[in] beta_z Argument in z-direction.
+ * \return Value of the integral.
+ */
+inline double GetBesselIntegral(double beta_x, double beta_y, double beta_z) {
+  const double A = 1.0 + 2.0 * (beta_x + beta_y + beta_z);
+  const double Bx = 2.0 * beta_x;
+  const double By = 2.0 * beta_y;
+  const double Bz = 2.0 * beta_z;
+  const int N = 4000;
+  const double t_max = 20.0;
+  const double dt = t_max / N;
+  double sum = 0.0;
+  for (int i = 1; i <= N; i++) {
+    double t = i * dt;
+    double lx = GetBesselZero(Bx * t);
+    double ly = GetBesselZero(By * t);
+    double lz = GetBesselZero(Bz * t);
+    double lin = log(t) - A * t + lx + ly + lz;
+    double integrand = exp(lin);
+    if (i==N) integrand *= 0.5;
+    sum += integrand;
+  }
+  return sum * dt;
+}
+
+/// @}
+}  // namespace RandomToolbox
\ No newline at end of file
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 8a95237ac3fb..3b1304c73e66 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -2960,9 +2960,48 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: DES Constant */
   addDoubleOption("DES_CONST", Const_DES, 0.65);
 
+  /* DESCRIPTION: SBS lengthscale coefficient */
+  addDoubleOption("SBS_LENGTHSCALE_COEFF", SBS_Cdelta, 0.02);
+
+  /* DESCRIPTION: Maximum number of smoothing iterations for SBS model. */
+  addUnsignedShortOption("SBS_MAX_ITER_SMOOTH", SBS_maxIterSmooth, 100);
+
+  /* DESCRIPTION: SBS timescale coefficient */
+  addDoubleOption("SBS_TIMESCALE_COEFF", SBS_Ctau, 0.05);
+
+  /* DESCRIPTION: SBS intensity coefficient */
+  addDoubleOption("SBS_INTENSITY_COEFF", SBS_Cmag, 1.0);
+
   /* DESCRIPTION: Specify Hybrid RANS/LES model */
   addEnumOption("HYBRID_RANSLES", Kind_HybridRANSLES, HybridRANSLES_Map, NO_HYBRIDRANSLES);
 
+  /* DESCRIPTION: Specify if the Stochastic Backscatter Model must be activated */
+  addBoolOption("STOCHASTIC_BACKSCATTER", StochasticBackscatter, false);
+
+  /* DESCRIPTION: Specify if the LES mode must be enforced */
+  addBoolOption("ENFORCE_LES", enforceLES, false);
+
+  /* DESCRIPTION: Specify if the stochastic source term must be included in the turbulence model equation */
+  addBoolOption("SBS_SOURCE_NU_EQUATION", stochSourceNu, false);
+
+  /* DESCRIPTION: Enable diagnostics of the stochastic source term in Langevin equation. */
+  addBoolOption("SBS_SOURCE_DIAGNOSTICS", stochSourceDiagnostics, false);
+
+  /* DESCRIPTION: Relaxation factor for the stochastic source term (Stochastic Backscatter Model) */
+  addDoubleOption("SBS_RELAXATION_FACTOR", stochSourceRelax, 0.0);
+
+  /* DESCRIPTION: Apply Stochastic Backscatter Model only in a bounded box */
+  addBoolOption("SBS_IN_BOX", StochBackscatterInBox, false);
+
+  /* DESCRIPTION: Specify extents of box where Stochastic Backscatter Model is active */
+  addDoubleArrayOption("SBS_BOX_BOUNDS", 6, false, StochBackscatterBoxBounds);
+
+  /* DESCRIPTION: Shielding function lower threshold for application of Stochastic Backscatter Model */
+  addDoubleOption("SBS_FD_LOWER_THRESHOLD", stochFdThreshold, 0.9);
+
+  /* DESCRIPTION: Filter width for LES (if negative, it is computed based on the local cell size) */
+  addDoubleOption("LES_FILTER_WIDTH", LES_FilterWidth, -1.0);
+
   /* DESCRIPTION: Roe with low dissipation for unsteady flows */
   addEnumOption("ROE_LOW_DISSIPATION", Kind_RoeLowDiss, RoeLowDiss_Map, NO_ROELOWDISS);
 
@@ -6522,6 +6561,58 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
           case SA_ZDES:  cout << "Delayed Detached Eddy Simulation (DDES) with Vorticity-based SGS" << endl; break;
           case SA_EDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break;
         }
+        if (Kind_HybridRANSLES != NO_HYBRIDRANSLES) {
+          if (LES_FilterWidth > 0.0) cout << "User-specified LES filter width: " << LES_FilterWidth << endl;
+          cout << "Stochastic Backscatter: ";
+          if (StochasticBackscatter) {
+            cout << "ON" << endl;
+            cout << "Backscatter intensity coefficient: " << SBS_Cmag << endl;
+            if (SBS_Cmag < 0.0)
+              SU2_MPI::Error("Backscatter intensity coefficient must be non-negative.", CURRENT_FUNCTION);
+            cout << "Backscatter lengthscale coefficient: " << SBS_Cdelta << endl;
+            if (SBS_Cdelta < 0.0)
+              SU2_MPI::Error("Backscatter lengthscale coefficient must be non-negative.", CURRENT_FUNCTION);
+            cout << "Backscatter timescale coefficient: " << SBS_Ctau << endl;
+            if (SBS_Ctau < 0.0)
+              SU2_MPI::Error("Backscatter timescale coefficient must be non-negative.", CURRENT_FUNCTION);
+            if (SBS_maxIterSmooth > 0)
+              cout << "Maximum number of iterations for implicit smoothing: " << SBS_maxIterSmooth << endl;
+            else
+              cout << "No smoothing applied to source terms in Langevin equations." << endl;
+            if (stochSourceNu)
+              cout << "Stochastic source term included in turbulence model equation." << endl;
+            else
+              cout << "Stochastic source term NOT included in turbulence model equation." << endl;
+            if (stochSourceRelax > 0.0)
+              cout << "Relaxation factor for stochastic source term: " << stochSourceRelax << endl;
+            else
+              cout << "No relaxation factor for stochastic source term." << endl;
+            if (StochBackscatterInBox) {
+              cout << "Stochastic Backscatter Model activated only in a bounded box." << endl;
+              cout << "Box bounds: " << endl;
+              cout << "  X: " << setw(10) << fixed << setprecision(4) << StochBackscatterBoxBounds[0] << " , "
+                              << setw(10) << fixed << setprecision(4) << StochBackscatterBoxBounds[1] << endl;
+              cout << "  Y: " << setw(10) << fixed << setprecision(4) << StochBackscatterBoxBounds[2] << " , "
+                              << setw(10) << fixed << setprecision(4) << StochBackscatterBoxBounds[3] << endl;
+              cout << "  Z: " << setw(10) << fixed << setprecision(4) << StochBackscatterBoxBounds[4] << " , "
+                              << setw(10) << fixed << setprecision(4) << StochBackscatterBoxBounds[5] << endl;
+            } else {
+              cout << "Stochastic Backscatter Model activated in the whole computational domain." << endl;
+            }
+            if (Kind_HybridRANSLES != SA_DES)
+              cout << "Stochastic source terms suppressed where the shielding function is lower than: " << setw(5) << setprecision(3) << stochFdThreshold << endl;
+          } else {
+            cout << "OFF" << endl;
+          }
+        }
+        if (StochasticBackscatter && Kind_HybridRANSLES == NO_HYBRIDRANSLES)
+          SU2_MPI::Error("Stochastic Backscatter can only be activated with Hybrid RANS/LES.", CURRENT_FUNCTION);
+        if (enforceLES) {
+          if (Kind_HybridRANSLES == NO_HYBRIDRANSLES)
+            SU2_MPI::Error("ENFORCE_LES can only be activated with Hybrid RANS/LES.", CURRENT_FUNCTION);
+          else
+            cout << "LES enforced in the whole computational domain." << endl;
+        }
         break;
       case MAIN_SOLVER::NEMO_EULER:
         if (Kind_Regime == ENUM_REGIME::COMPRESSIBLE) cout << "Compressible two-temperature thermochemical non-equilibrium Euler equations." << endl;
@@ -7090,11 +7181,20 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
           cout << "Lax viscous coefficients (1st): " << Kappa_1st_Flow << ".\n";
           cout << "First order integration." << endl;
         }
+        else if (Kind_Centered_Flow == CENTERED::LD2) {
+          cout << "Low-Dissipation Low-Dispersion (LD2) scheme for the flow inviscid terms." << endl;
+          if (!(Kind_Solver==MAIN_SOLVER::INC_EULER || Kind_Solver==MAIN_SOLVER::INC_NAVIER_STOKES || Kind_Solver==MAIN_SOLVER::INC_RANS))
+            SU2_MPI::Error("LD2 scheme not yet implemented for the compressible flow solver.", CURRENT_FUNCTION);
+          if (Kind_FluidModel != CONSTANT_DENSITY)
+            SU2_MPI::Error("LD2 scheme available for constant density flows only.", CURRENT_FUNCTION);
+          if (Energy_Equation)
+            cout << "WARNING: Current implementation of the LD2 scheme not compatible with the energy equation. JST employed in energy equation instead." << endl;
+        }
         else {
-          cout << "Jameson-Schmidt-Turkel scheme (2nd order in space) for the flow inviscid terms.\n";
-          cout << "JST viscous coefficients (2nd & 4th): " << Kappa_2nd_Flow << ", " << Kappa_4th_Flow << ".\n";
-          cout << "The method includes a grid stretching correction (p = 0.3)."<< endl;
+            cout << "Jameson-Schmidt-Turkel scheme (2nd order in space) for the flow inviscid terms.\n";
         }
+        cout << "JST viscous coefficients (2nd & 4th): " << Kappa_2nd_Flow << ", " << Kappa_4th_Flow << ".\n";
+        cout << "The method includes a grid stretching correction (p = 0.3)."<< endl;
       }
 
       if (Kind_ConvNumScheme_Flow == SPACE_UPWIND) {
diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp
index cc9194997478..d5c11569698c 100644
--- a/SU2_CFD/include/numerics/CNumerics.hpp
+++ b/SU2_CFD/include/numerics/CNumerics.hpp
@@ -32,6 +32,7 @@
 #include <iostream>
 #include <limits>
 #include <cstdlib>
+#include <random>
 
 #include "../../../Common/include/CConfig.hpp"
 #include "../../../Common/include/linear_algebra/blas_structure.hpp"
@@ -184,6 +185,14 @@ class CNumerics {
   roughness_j = 0.0;                       /*!< \brief Roughness of the wall nearest to point j. */
 
   su2double MeanPerturbedRSM[3][3];   /*!< \brief Perturbed Reynolds stress tensor  */
+  su2double stochReynStress[3][3] = {{0.0}}; /*!< \brief Stochastic contribution to Reynolds stress tensor for Backscatter Model. */
+  su2double stochSource[3] = {0.0}; /*!< \brief Source term for Langevin equations in Stochastic Backscatter Model. */
+  su2double
+  stochVar_i[3] = {0.0}, /*!< \brief Stochastic variables at point i for Stochastic Backscatter Model. */
+  stochVar_j[3] = {0.0}; /*!< \brief Stochastic variables at point j for Stochastic Backscatter Model. */
+  su2double
+  lesMode_i = 0.0, /*!< \brief LES sensor at point i for hybrid RANS-LES methods. */
+  lesMode_j = 0.0; /*!< \brief LES sensor at point j for hybrid RANS-LES methods. */
   SST_ParsedOptions sstParsedOptions; /*!< \brief additional options for the SST turbulence model */
   unsigned short Eig_Val_Comp;    /*!< \brief Component towards which perturbation is perfromed */
   su2double uq_delta_b;           /*!< \brief Magnitude of perturbation */
@@ -640,6 +649,63 @@ class CNumerics {
       }
     }
   }
+  
+  /*!
+   * \brief Compute a random contribution to the Reynolds stress tensor (Stochastic Backscatter Model).
+   * \details See: Kok, Johan C. "A stochastic backscatter model for grey-area mitigation in detached
+   * eddy simulations." Flow, Turbulence and Combustion 99.1 (2017): 119-150.
+   * \param[in] nDim - Dimension of the flow problem, 2 or 3.
+   * \param[in] density - Density.
+   * \param[in] tke - Turbulent kinetic energy.
+   * \param[in] rndVec - Vector of stochastic variables from Langevin equations.
+   * \param[in] Cmag - Stochastic backscatter intensity coefficient.
+   * \param[out] stochReynStress - Stochastic tensor (to be added to the Reynolds stress tensor).
+   */
+  template<class Mat, class Scalar, class Vector>
+  NEVERINLINE static void ComputeStochReynStress(size_t nDim, Scalar density, Scalar tke,
+                                                 Vector rndVec, Scalar Cmag, Mat& stochReynStress) {
+
+    /* --- Calculate stochastic tensor --- */
+
+    su2double stochLim = 3.0; 
+
+    stochReynStress[0][0] = 0.0;
+    stochReynStress[1][1] = 0.0;
+    stochReynStress[2][2] = 0.0;
+    stochReynStress[0][1] =   Cmag * density * tke * max(-stochLim, min(stochLim, rndVec[2]));
+    stochReynStress[0][2] = - Cmag * density * tke * max(-stochLim, min(stochLim, rndVec[1]));
+    stochReynStress[1][2] =   Cmag * density * tke * max(-stochLim, min(stochLim, rndVec[0]));
+    stochReynStress[1][0] = - stochReynStress[1][0];
+    stochReynStress[2][0] = - stochReynStress[2][0];
+    stochReynStress[2][1] = - stochReynStress[2][1];
+
+  }
+
+    /*!
+   * \brief Compute relaxation factor for stochastic source term in momentum equations (Stochastic Backscatter Model).
+   * \param[in] config - Definition of the particular problem.
+   * \param[out] intensityCoeff - Relaxation factor for backscatter intensity.
+   */
+  NEVERINLINE static su2double ComputeStochRelaxFactor(const CConfig* config) {
+
+    unsigned long timeIter = config->GetTimeIter();
+    unsigned long restartIter = config->GetRestart_Iter();
+    su2double SBS_Cmag = config->GetSBS_Cmag();
+    su2double SBS_RelaxFactor = config->GetStochSourceRelax();
+    su2double intensityCoeff = SBS_Cmag;
+    if (SBS_RelaxFactor > 0.0) {
+      su2double FS_Vel = config->GetModVel_FreeStream();
+      su2double ReynoldsLength = config->GetLength_Reynolds();
+      su2double timeScale = ReynoldsLength / FS_Vel;
+      unsigned long timeIter = config->GetTimeIter();
+      unsigned long restartIter = config->GetRestart_Iter();
+      su2double timeStep = config->GetTime_Step();
+      su2double currentTime = (timeIter - restartIter) * timeStep;
+      intensityCoeff = SBS_Cmag * (1.0 - exp(- currentTime / (timeScale*SBS_RelaxFactor)));
+    }
+    return intensityCoeff;
+
+  }
 
   /*!
    * \brief Project average gradient onto normal (with or w/o correction) for viscous fluxes of scalar quantities.
@@ -838,6 +904,27 @@ class CNumerics {
     turb_ke_j = val_turb_ke_j;
   }
 
+  /*!
+   * \brief Set the stochastic variables from Langevin equations (Stochastic Backscatter Model).
+   * \param[in] val_stochvar_i - Value of the stochastic variable at point i.
+   * \param[in] val_stochvar_j - Value of the stochastic variable at point j.
+   * \param[in] iDim - Index of Langevin equation.
+   */
+  inline void SetStochVar(unsigned short iDim, su2double val_stochvar_i, su2double val_stochvar_j) {
+    stochVar_i[iDim] = val_stochvar_i;
+    stochVar_j[iDim] = val_stochvar_j;
+  }
+
+  /*!
+   * \brief Set the LES sensor for hybrid RANS-LES methods.
+   * \param[in] val_lesMode_i - Value of the LES sensor at point i.
+   * \param[in] val_lesMode_j - Value of the LES sensor at point j.
+   */
+  inline void SetLES_Mode(su2double val_lesMode_i, su2double val_lesMode_j) {
+    lesMode_i = val_lesMode_i;
+    lesMode_j = val_lesMode_j;
+  }
+
   /*!
    * \brief Set the value of the distance from the nearest wall.
    * \param[in] val_dist_i - Value of of the distance from point i to the nearest wall.
@@ -848,6 +935,15 @@ class CNumerics {
     dist_j = val_dist_j;
   }
 
+  /*!
+   * \brief Set the stochastic source term for the Langevin equations (Backscatter Model).
+   * \param[in] val_stoch_source - Value of stochastic source term.
+   * \param[in] iDim - Index of Langevin equation.
+   */
+  void SetStochSource(su2double val_stoch_source, unsigned short iDim) {
+    stochSource[iDim] = val_stoch_source;
+  }
+
   /*!
    * \brief Set the value of the roughness from the nearest wall.
    * \param[in] val_dist_i - Value of of the roughness of the nearest wall from point i
diff --git a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp
index 01d5342abb42..57dfea12ff25 100644
--- a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp
+++ b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp
@@ -61,6 +61,7 @@ class CAvgGrad_Base : public CNumerics {
   Mean_Cp,                                /*!< \brief Mean value of the specific heat capacity at constant pressure. */
   Mean_turb_ke,                           /*!< \brief Mean value of the turbulent kinetic energy. */
   Mean_TauWall,                           /*!< \brief Mean wall shear stress (wall functions). */
+  Mean_StochVar[3],                       /*!< \brief Mean stochastic variables (Stochastic Backscatter Model). */
   TauWall_i, TauWall_j,                   /*!< \brief Wall shear stress at point i and j (wall functions). */
   dist_ij_2,                              /*!< \brief Length of the edge and face, squared */
   Edge_Vector[MAXNDIM] = {0.0};           /*!< \brief Vector from point i to point j. */
@@ -199,12 +200,14 @@ class CAvgGrad_Base : public CNumerics {
    * \param[in] val_turb_ke - Turbulent kinetic energy
    * \param[in] val_laminar_viscosity - Laminar viscosity.
    * \param[in] val_eddy_viscosity - Eddy viscosity.
+   * \param[in] config - Definition of the particular problem.
    */
   void SetStressTensor(const su2double *val_primvar,
                        const su2double* const *val_gradprimvar,
                        su2double val_turb_ke,
                        su2double val_laminar_viscosity,
-                       su2double val_eddy_viscosity);
+                       su2double val_eddy_viscosity,
+                       const CConfig* config);
 
   /*!
    * \brief Get a component of the viscous stress tensor.
diff --git a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp
index b8d8a6403930..41d0ca2c5ac9 100644
--- a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp
+++ b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp
@@ -52,6 +52,7 @@ class CUpwScalar : public CNumerics {
   const FlowIndices idx;            /*!< \brief Object to manage the access to the flow primitives. */
   su2double a0 = 0.0;               /*!< \brief The maximum of the face-normal velocity and 0. */
   su2double a1 = 0.0;               /*!< \brief The minimum of the face-normal velocity and 0. */
+  su2double m_ij = 0.0;              /*!< \brief Face-normal momentum (Langevin equations). */
   su2double Flux[MAXNVAR];          /*!< \brief Final result, diffusive flux/residual. */
   su2double* Jacobian_i[MAXNVAR];   /*!< \brief Flux Jacobian w.r.t. node i. */
   su2double* Jacobian_j[MAXNVAR];   /*!< \brief Flux Jacobian w.r.t. node j. */
@@ -140,6 +141,13 @@ class CUpwScalar : public CNumerics {
       a1 = fmin(0.0, q_ij);
     }
 
+    if (config->GetStochastic_Backscatter()) {
+      m_ij = 0.0;
+      for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+        m_ij += 0.5 * (V_i[idx.Density()]*V_i[iDim + idx.Velocity()] + V_j[idx.Density()]*V_j[iDim + idx.Velocity()]) * Normal[iDim];
+      }
+    }
+
     FinishResidualCalc(config);
 
     AD::SetPreaccOut(Flux, nVar);
diff --git a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp
index a2942fc72999..44256ad04058 100644
--- a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp
+++ b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp
@@ -64,7 +64,7 @@ class CAvgGrad_Scalar : public CNumerics {
   const FlowIndices idx;                      /*!< \brief Object to manage the access to the flow primitives. */
   su2double Proj_Mean_GradScalarVar[MAXNVAR]; /*!< \brief Mean_gradScalarVar DOT normal, corrected if required. */
   su2double proj_vector_ij = 0.0;             /*!< \brief (Edge_Vector DOT normal)/|Edge_Vector|^2 */
-  su2double Flux[MAXNVAR];                    /*!< \brief Final result, diffusive flux/residual. */
+  su2double Flux[MAXNVAR] = {0.0};            /*!< \brief Final result, diffusive flux/residual. */
   su2double* Jacobian_i[MAXNVAR];             /*!< \brief Flux Jacobian w.r.t. node i. */
   su2double* Jacobian_j[MAXNVAR];             /*!< \brief Flux Jacobian w.r.t. node j. */
   su2double JacobianBuffer[2*MAXNVAR*MAXNVAR];/*!< \brief Static storage for the two Jacobians. */
diff --git a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp
index 5e401210734a..2df13c7bd55c 100644
--- a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp
+++ b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp
@@ -48,6 +48,11 @@ class CUpwSca_TurbSA final : public CUpwScalar<FlowIndices> {
   using Base::ScalarVar_i;
   using Base::ScalarVar_j;
   using Base::bounded_scalar;
+  using Base::V_i;
+  using Base::V_j;
+  using Base::idx;
+  using Base::nVar;
+  using Base::m_ij;
 
   /*!
    * \brief Adds any extra variables to AD.
@@ -59,6 +64,17 @@ class CUpwSca_TurbSA final : public CUpwScalar<FlowIndices> {
    * \param[in] config - Definition of the particular problem.
    */
   void FinishResidualCalc(const CConfig* config) override {
+    if (config->GetStochastic_Backscatter()) {
+      for (unsigned short iVar = 1; iVar < nVar; iVar++) {
+        Flux[iVar] = m_ij * 0.5 * (ScalarVar_i[iVar] + ScalarVar_j[iVar]);
+      }
+      for (unsigned short iVar = 0; iVar < nVar; iVar++) {
+        for (unsigned short jVar = 0; jVar < nVar; jVar++) {
+          Jacobian_i[iVar][jVar] = (iVar == jVar) ? 0.5*m_ij : 0.0;
+          Jacobian_j[iVar][jVar] = (iVar == jVar) ? 0.5*m_ij : 0.0;
+        }
+      }
+    }
     Flux[0] = a0*ScalarVar_i[0] + a1*ScalarVar_j[0];
     Jacobian_i[0][0] = a0;
     Jacobian_j[0][0] = a1;
diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
index a5217a081df3..e9462a572741 100644
--- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
+++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
@@ -72,8 +72,8 @@ class CSourceBase_TurbSA : public CNumerics {
  protected:
 
   /*--- Residual and Jacobian ---*/
-  su2double Residual, *Jacobian_i;
-  su2double Jacobian_Buffer; /*!< \brief Static storage for the Jacobian (which needs to be pointer for return type). */
+  su2double Residual[4], *Jacobian_i[4]; /*!< \brief Increase the size of residual and Jacobian for Langevin equations (Stochastic Backscatter Model).*/
+  su2double Jacobian_Buffer[16]; /*!< \brief Static storage for the Jacobian (which needs to be pointer for return type). */ 
 
   const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */
   const SA_ParsedOptions options; /*!< \brief Struct with SA options. */
@@ -109,7 +109,79 @@ class CSourceBase_TurbSA : public CNumerics {
     /* Diffusion source term */
     const su2double dv_axi = (1.0/sigma)*nu_e*ScalarVar_Grad_i[0][1];
 
-    Residual += yinv * dv_axi * Volume;
+    Residual[0] += yinv * dv_axi * Volume;
+  }
+
+  /*!
+   * \brief Include source-term residuals for Langevin equations (Stochastic Backscatter Model) 
+   */
+  inline void ResidualStochEquations(su2double timeStep, const su2double ct, 
+                                     su2double lengthScale, su2double DES_const,
+                                     const CSAVariables& var, TIME_MARCHING time_marching,
+                                     su2double threshold) {
+
+    const su2double& nue = ScalarVar_i[0];
+    const auto& density = V_i[idx.Density()];
+    const su2double limiterRANS = 0.1;
+    const su2double timeStepFrac = 0.2;
+    const su2double nut = max(nue*var.fv1, 1e-10);
+    const su2double delta = lengthScale/DES_const;
+
+    if (delta > 1e-10) {
+
+      su2double tTurb = ct*pow(delta, 2)/nut;
+      su2double blendingFactor = DES_const*DES_const/ct*(1.0-lesMode_i) + lesMode_i;
+      tTurb *= blendingFactor;
+      tTurb = max(timeStep*timeStepFrac, min(tTurb, timeStep/timeStepFrac));
+      su2double tRat = timeStep / tTurb;
+    
+      su2double corrFac = 1.0;
+      if (time_marching == TIME_MARCHING::DT_STEPPING_2ND) {
+        corrFac = sqrt(0.5*(1.0+tRat)*(4.0+tRat)/(2.0+tRat));
+      } else if (time_marching == TIME_MARCHING::DT_STEPPING_1ST) {
+        corrFac = sqrt(1.0+0.5*tRat);
+      }
+    
+      su2double scaleFactor = 0.0;
+      if (lesMode_i > threshold) scaleFactor = 1.0/tTurb * sqrt(2.0/tRat) * density * corrFac;
+
+      for (unsigned short iVar = 1; iVar < nVar; iVar++) {
+        Residual[iVar] = scaleFactor * stochSource[iVar-1] - 1.0/tTurb * density * ScalarVar_i[iVar];
+        Residual[iVar] *= Volume;
+      }
+
+      for (unsigned short iVar = 1; iVar < nVar; iVar++ )
+        Jacobian_i[iVar][iVar] = -1.0/tTurb * density * Volume;
+
+    }
+
+  }
+
+  /*!
+   * \brief Include stochastic source term in the Spalart-Allmaras turbulence model equation (Stochastic Backscatter Model).
+   */
+  inline void AddStochSource(CSAVariables& var, const su2double Cmag, su2double threshold, su2double& prod) {
+
+    su2double nut = ScalarVar_i[0] * var.fv1;
+    su2double tke = 0.0;
+    const su2double limiter = 5.0;
+    if (lesMode_i > threshold) tke = nut * StrainMag_i;
+
+    su2double R12 =   Cmag * tke * ScalarVar_i[3];
+    su2double R13 = - Cmag * tke * ScalarVar_i[2];
+    su2double R23 =   Cmag * tke * ScalarVar_i[1];
+
+    su2double RGradU = R12*Vorticity_i[2] - R13*Vorticity_i[1] + R23*Vorticity_i[0];
+
+    su2double stochProdK = - RGradU;
+    su2double Ji_3 = pow(var.Ji, 3);
+    su2double Dfv1Dnut = 3.0 * var.fv1 * var.cv1_3 / (var.cv1_3 + Ji_3);
+    su2double fac = 1.0 / (var.fv1 + ScalarVar_i[0]*Dfv1Dnut);
+    su2double stochProdNut = stochProdK * dist_i*dist_i/(2.0*ScalarVar_i[0]) * fac;
+    stochProdNut = max(-limiter*prod, min(limiter*prod, stochProdNut));
+
+    prod += stochProdNut;
+
   }
 
  public:
@@ -119,17 +191,19 @@ class CSourceBase_TurbSA : public CNumerics {
    * \param[in] config - Definition of the particular problem.
    */
   CSourceBase_TurbSA(unsigned short nDim, const CConfig* config)
-      : CNumerics(nDim, 1, config),
+      : CNumerics(nDim, config->GetStochastic_Backscatter() ? 4 : 1, config),
         idx(nDim, config->GetnSpecies()),
         options(config->GetSAParsedOptions()),
         axisymmetric(config->GetAxisymmetric()),
         transition_LM(config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) {
     /*--- Setup the Jacobian pointer, we need to return su2double** but we know
      * the Jacobian is 1x1 so we use this trick to avoid heap allocation. ---*/
-    Jacobian_i = &Jacobian_Buffer;
+    //Jacobian_i = &Jacobian_Buffer;
+    /*--- Setup the Jacobian pointer (size increased for Stochastic Backscatter Model). ---*/
+    for (unsigned short iVar = 0; iVar < 4; iVar++)
+      Jacobian_i[iVar] = Jacobian_Buffer + 4*iVar;
   }
 
-
   /*!
    * \brief Residual for source term integration.
    * \param[in] config - Definition of the particular problem.
@@ -140,16 +214,22 @@ class CSourceBase_TurbSA : public CNumerics {
     const auto& laminar_viscosity = V_i[idx.LaminarViscosity()];
 
     AD::StartPreacc();
-    AD::SetPreaccIn(density, laminar_viscosity, StrainMag_i, ScalarVar_i[0], Volume, dist_i, roughness_i);
+    AD::SetPreaccIn(density, laminar_viscosity, StrainMag_i, Volume, dist_i, roughness_i);
+    AD::SetPreaccIn(ScalarVar_i, nVar);
     AD::SetPreaccIn(Vorticity_i, 3);
     AD::SetPreaccIn(PrimVar_Grad_i + idx.Velocity(), nDim, nDim);
-    AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim);
+    AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim);
+    AD::SetPreaccIn(stochSource, 3);
 
     /*--- Common auxiliary variables and constants of the model. ---*/
     CSAVariables var;
 
-    Residual = 0.0;
-    Jacobian_i[0] = 0.0;
+    for (unsigned short iVar = 0; iVar < 4; iVar++) {
+      Residual[iVar] = 0.0;
+      for (unsigned short jVar = 0; jVar < 4; jVar++) {
+        Jacobian_i[iVar][jVar] = 0.0;
+      }
+    }
 
     if (dist_i > 1e-10) {
 
@@ -245,19 +325,33 @@ class CSourceBase_TurbSA : public CNumerics {
 
       /*--- Compute production, destruction and jacobian ---*/
       su2double Production = 0.0, Destruction = 0.0;
-      SourceTerms::get(ScalarVar_i[0], var, Production, Destruction, Jacobian_i[0]);
+      SourceTerms::get(ScalarVar_i[0], var, Production, Destruction, Jacobian_i[0][0]);
+
+      if (config->GetStochastic_Backscatter() && config->GetStochSourceNu()) {
+        su2double intensityCoeff = ComputeStochRelaxFactor(config);
+        AddStochSource(var, intensityCoeff, config->GetStochFdThreshold(), Production);
+      }
 
-      Residual = (Production - Destruction) * Volume;
+      Residual[0] = (Production - Destruction) * Volume;
 
       if (axisymmetric) ResidualAxisymmetricDiffusion(var.sigma);
 
-      Jacobian_i[0] *= Volume;
+      Jacobian_i[0][0] *= Volume;
+
+      /*--- Compute residual for Langevin equations (Stochastic Backscatter Model). ---*/
+
+      if (config->GetStochastic_Backscatter()) {
+        const su2double DES_const = config->GetConst_DES();
+        const su2double ctau = config->GetSBS_Ctau();
+        ResidualStochEquations(config->GetDelta_UnstTime(), ctau, dist_i, DES_const,
+                               var, config->GetTime_Marching(), config->GetStochFdThreshold());
+      }
     }
 
-    AD::SetPreaccOut(Residual);
+    AD::SetPreaccOut(Residual, 4);
     AD::EndPreacc();
 
-    return ResidualType<>(&Residual, &Jacobian_i, nullptr);
+    return ResidualType<>(Residual, Jacobian_i, nullptr);
   }
 };
 
@@ -553,14 +647,14 @@ class CCompressibilityCorrection final : public ParentClass {
       const su2double d_axiCorrection = 2.0 * c5 * nue * pow(v * yinv / sound_speed, 2) * Volume;
       const su2double axiCorrection = 0.5 * nue * d_axiCorrection;
 
-      this->Residual -= axiCorrection;
-      this->Jacobian_i[0] -= d_axiCorrection;
+      this->Residual[0] -= axiCorrection;
+      this->Jacobian_i[0][0] -= d_axiCorrection;
     }
 
-    this->Residual -= CompCorrection;
-    this->Jacobian_i[0] -= d_CompCorrection;
+    this->Residual[0] -= CompCorrection;
+    this->Jacobian_i[0][0] -= d_CompCorrection;
 
-    return ResidualType(&this->Residual, &this->Jacobian_i, nullptr);
+    return ResidualType(this->Residual, this->Jacobian_i, nullptr);
   }
 };
 
diff --git a/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp b/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp
index 09e71002f8ae..69101654c99c 100644
--- a/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp
+++ b/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp
@@ -71,6 +71,7 @@ CNumericsSIMD* createCenteredNumerics(const CConfig& config, int iMesh, const CV
       obj = new CLaxScheme<ViscousDecorator>(config, iMesh, turbVars);
       break;
     case CENTERED::JST:
+    case CENTERED::LD2: // Just to silence compiler warnings (LD2 implemented in the incompressible solver only).
       obj = new CJSTScheme<ViscousDecorator>(config, iMesh, turbVars);
       break;
     case CENTERED::JST_KE:
diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
index 8e935770a89d..5129e9b8b8ec 100644
--- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
+++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
@@ -434,9 +434,10 @@ void CFVMFlowSolverBase<V, R>::Viscous_Residual_impl(unsigned long iEdge, CGeome
 
   const bool implicit  = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
   const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST);
+  const bool backscatter = config->GetStochastic_Backscatter();
 
   CVariable* turbNodes = nullptr;
-  if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes();
+  if (tkeNeeded || backscatter) turbNodes = solver_container[TURB_SOL]->GetNodes();
 
   /*--- Points, coordinates and normal vector in edge ---*/
 
@@ -467,6 +468,28 @@ void CFVMFlowSolverBase<V, R>::Viscous_Residual_impl(unsigned long iEdge, CGeome
     numerics->SetTurbKineticEnergy(turbNodes->GetSolution(iPoint,0),
                                    turbNodes->GetSolution(jPoint,0));
 
+  /*--- Stochastic variables from Langevin equations (Stochastic Backscatter Model). ---*/
+
+  if (backscatter) {
+    for (unsigned short iDim = 0; iDim < nDim; iDim++)
+      numerics->SetStochVar(iDim, turbNodes->GetSolution(iPoint, iDim+1),
+                                  turbNodes->GetSolution(jPoint, iDim+1));
+    su2double DES_length_i = turbNodes->GetDES_LengthScale(iPoint);
+    su2double DES_length_j = turbNodes->GetDES_LengthScale(jPoint);
+    su2double lesMode_i = (DES_length_i > 1e-10) ? turbNodes->GetLES_Mode(iPoint) : 0.0;
+    su2double lesMode_j = (DES_length_j > 1e-10) ? turbNodes->GetLES_Mode(jPoint) : 0.0;
+    su2double tke_i = 0.0, tke_j = 0.0;
+    if (max(lesMode_i, lesMode_j) > config->GetStochFdThreshold()) {
+      su2double eddyVisc_i = turbNodes->GetmuT(iPoint) / nodes->GetDensity(iPoint);
+      su2double eddyVisc_j = turbNodes->GetmuT(jPoint) / nodes->GetDensity(jPoint);
+      su2double strainMag_i = nodes->GetStrainMag(iPoint);
+      su2double strainMag_j = nodes->GetStrainMag(jPoint);
+      tke_i = strainMag_i * eddyVisc_i;
+      tke_j = strainMag_j * eddyVisc_j;
+    }
+    numerics->SetTurbKineticEnergy(tke_i, tke_j);
+  }
+
   /*--- Wall shear stress values (wall functions) ---*/
 
   numerics->SetTau_Wall(nodes->GetTau_Wall(iPoint),
diff --git a/SU2_CFD/include/solvers/CTurbSASolver.hpp b/SU2_CFD/include/solvers/CTurbSASolver.hpp
index dfc1362b003c..a19ffb3dd9e4 100644
--- a/SU2_CFD/include/solvers/CTurbSASolver.hpp
+++ b/SU2_CFD/include/solvers/CTurbSASolver.hpp
@@ -39,7 +39,9 @@
 
 class CTurbSASolver final : public CTurbSolver {
 private:
-  su2double nu_tilde_Engine, nu_tilde_ActDisk;
+
+  su2double nu_tilde_Engine[4] = {0.0};
+  su2double nu_tilde_ActDisk[4] = {0.0};
 
   /*!
    * \brief A virtual member.
@@ -51,6 +53,20 @@ class CTurbSASolver final : public CTurbSolver {
                           CGeometry *geometry,
                           CConfig *config);
 
+  /*!
+   * \brief Update the source terms of the stochastic equations (Stochastic Backscatter Model).
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] geometry - Geometrical definition.
+   */
+  void SetLangevinSourceTerms(CConfig *config, CGeometry* geometry);
+
+  /*!
+   * \brief Apply Laplacian smoothing to the source terms in Langevin equations (Stochastic Backscatter Model).
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] geometry - Geometrical definition.
+   */
+  void SmoothLangevinSourceTerms(CConfig* config, CGeometry* geometry);
+
   /*!
    * \brief Compute nu tilde from the wall functions.
    * \param[in] geometry - Geometrical definition of the problem.
diff --git a/SU2_CFD/include/variables/CIncNSVariable.hpp b/SU2_CFD/include/variables/CIncNSVariable.hpp
index 7dd082cdca28..5aa1ebaa75e9 100644
--- a/SU2_CFD/include/variables/CIncNSVariable.hpp
+++ b/SU2_CFD/include/variables/CIncNSVariable.hpp
@@ -39,7 +39,8 @@
 class CIncNSVariable final : public CIncEulerVariable {
 private:
   VectorType Tau_Wall;        /*!< \brief Magnitude of the wall shear stress from a wall function. */
-  VectorType DES_LengthScale;
+  VectorType DES_LengthScale; /*!< \brief DES Length Scale. */
+  VectorType lesMode;        /*!< \brief Sensor for local simulation mode (0=RANS, 1=LES).*/
   const bool Energy;          /*!< \brief Flag for Energy equation in incompressible flows. */
 
 public:
@@ -133,4 +134,17 @@ class CIncNSVariable final : public CIncEulerVariable {
    */
   inline su2double GetDES_LengthScale(unsigned long iPoint) const override { return DES_LengthScale(iPoint); }
 
+  /*!
+   * \brief Set the LES sensor.
+   */
+  inline void SetLES_Mode(unsigned long iPoint, su2double val_les_mode) override {
+    lesMode(iPoint) = val_les_mode;
+  }
+
+  /*!
+   * \brief Get the LES sensor.
+   * \return Value of the LES sensor.
+   */
+  inline su2double GetLES_Mode(unsigned long iPoint) const override { return lesMode(iPoint); }
+
 };
diff --git a/SU2_CFD/include/variables/CNEMONSVariable.hpp b/SU2_CFD/include/variables/CNEMONSVariable.hpp
index 9ec9e0fd00fd..4585c1c6b9b6 100644
--- a/SU2_CFD/include/variables/CNEMONSVariable.hpp
+++ b/SU2_CFD/include/variables/CNEMONSVariable.hpp
@@ -52,6 +52,7 @@ class CNEMONSVariable final : public CNEMOEulerVariable {
 
   VectorType Tau_Wall;          /*!< \brief Magnitude of the wall shear stress from a wall function. */
   VectorType DES_LengthScale;   /*!< \brief DES Length Scale. */
+  VectorType lesMode;          /*!< \brief Sensor for local simulation mode (0=RANS, 1=LES). */
   VectorType Roe_Dissipation;   /*!< \brief Roe low dissipation coefficient. */
   VectorType Vortex_Tilting;    /*!< \brief Value of the vortex tilting variable for DES length scale computation. */
 
diff --git a/SU2_CFD/include/variables/CNSVariable.hpp b/SU2_CFD/include/variables/CNSVariable.hpp
index ea8b263fd98d..9c3ab8dabf50 100644
--- a/SU2_CFD/include/variables/CNSVariable.hpp
+++ b/SU2_CFD/include/variables/CNSVariable.hpp
@@ -41,6 +41,7 @@ class CNSVariable final : public CEulerVariable {
 
   VectorType Tau_Wall;        /*!< \brief Magnitude of the wall shear stress from a wall function. */
   VectorType DES_LengthScale; /*!< \brief DES Length Scale. */
+  VectorType lesMode;        /*!< \brief Sensor for local simulation mode (0=RANS, 1=LES).*/
   VectorType Roe_Dissipation; /*!< \brief Roe low dissipation coefficient. */
   VectorType Vortex_Tilting;  /*!< \brief Value of the vortex tilting variable for DES length scale computation. */
 
@@ -185,6 +186,19 @@ class CNSVariable final : public CEulerVariable {
     DES_LengthScale(iPoint) = val_des_lengthscale;
   }
 
+/*!
+   * \brief Set the LES sensor.
+   */
+  inline void SetLES_Mode(unsigned long iPoint, su2double val_les_mode) override {
+    lesMode(iPoint) = val_les_mode;
+  }
+
+  /*!
+   * \brief Get the LES sensor.
+   * \return Value of the LES sensor.
+   */
+  inline su2double GetLES_Mode(unsigned long iPoint) const override { return lesMode(iPoint); }
+
   /*!
    * \brief Set the new solution for Roe Dissipation.
    * \param[in] val_delta - A scalar measure of the grid size
diff --git a/SU2_CFD/include/variables/CTurbSAVariable.hpp b/SU2_CFD/include/variables/CTurbSAVariable.hpp
index 3f55883ce9ca..7e4d88df4d2c 100644
--- a/SU2_CFD/include/variables/CTurbSAVariable.hpp
+++ b/SU2_CFD/include/variables/CTurbSAVariable.hpp
@@ -28,6 +28,7 @@
 #pragma once
 
 #include "CTurbVariable.hpp"
+#include <random>
 
 /*!
  * \class CTurbSAVariable
@@ -40,7 +41,11 @@ class CTurbSAVariable final : public CTurbVariable {
 
 private:
   VectorType DES_LengthScale;
+  VectorType lesMode;
+  MatrixType stochSource;
+  MatrixType stochSourceOld;
   VectorType Vortex_Tilting;
+  VectorType besselIntegral;
 
 public:
   /*!
@@ -73,6 +78,49 @@ class CTurbSAVariable final : public CTurbVariable {
    */
   inline void SetDES_LengthScale(unsigned long iPoint, su2double val_des_lengthscale) override { DES_LengthScale(iPoint) = val_des_lengthscale; }
 
+  /*!
+   * \brief Get the source terms for the stochastic equations.
+   * \param[in] iPoint - Point index.
+   * \param[in] iDim - Dimension index.
+   * \return Value of the source term for the stochastic equations.
+   */
+  inline su2double GetLangevinSourceTerms(unsigned long iPoint, unsigned short iDim) const override { return stochSource(iPoint, iDim); }
+
+  /*!
+   * \brief Set the source terms for the stochastic equations.
+   * \param[in] iPoint - Point index.
+   * \param[in] iDim - Dimension index.
+   * \param[in] val_stochSource - Value of the source term for the stochastic equations.
+   */
+  inline void SetLangevinSourceTerms(unsigned long iPoint, unsigned short iDim, su2double val_stochSource) override { stochSource(iPoint, iDim) = val_stochSource; }
+
+  /*!
+   * \brief Get the old value of the source terms for the stochastic equations.
+   * \param[in] iPoint - Point index.
+   * \param[in] iDim - Dimension index.
+   * \return Old value of the source terms for the stochastic equations.
+   */
+  inline su2double GetLangevinSourceTermsOld(unsigned long iPoint, unsigned short iDim) const override { return stochSourceOld(iPoint, iDim); }
+
+  /*!
+   * \brief Set the old value of source terms for the stochastic equations.
+   * \param[in] iPoint - Point index.
+   * \param[in] iDim - Dimension index.
+   * \param[in] val_stochSource_old - Old value of the source term for the stochastic equations.
+   */
+  inline void SetLangevinSourceTermsOld(unsigned long iPoint, unsigned short iDim, su2double val_stochSource_old) override { stochSourceOld(iPoint, iDim) = val_stochSource_old; }
+
+  /*!
+   * \brief Set the LES sensor.
+   */
+  inline void SetLES_Mode(unsigned long iPoint, su2double val_les_mode) override { lesMode(iPoint) = val_les_mode; }
+
+  /*!
+   * \brief Get the LES sensor.
+   * \return Value of the LES sensor.
+   */
+  inline su2double GetLES_Mode(unsigned long iPoint) const override { return lesMode(iPoint); }
+
   /*!
    * \brief Set the vortex tilting measure for computation of the EDDES length scale
    * \param[in] iPoint - Point index.
@@ -87,4 +135,17 @@ class CTurbSAVariable final : public CTurbVariable {
    */
   inline su2double GetVortex_Tilting(unsigned long iPoint) const override { return Vortex_Tilting(iPoint); }
 
+  /*!
+   * \brief Set the integral of the product of three Bessel functions appearing in Laplacian smoothing.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_integral - Value of the integral.
+   */
+  inline void SetBesselIntegral(unsigned long iPoint, su2double val_integral) override { besselIntegral(iPoint) = val_integral; }
+
+  /*!
+   * \brief Get the the integral of the product of three Bessel functions appearing in Laplacian smoothing.
+   * \return Value of the integral.
+   */
+  inline su2double GetBesselIntegral(unsigned long iPoint) const override { return besselIntegral(iPoint); }
+
 };
diff --git a/SU2_CFD/include/variables/CTurbVariable.hpp b/SU2_CFD/include/variables/CTurbVariable.hpp
index b4e4ae569743..34bf86e04e76 100644
--- a/SU2_CFD/include/variables/CTurbVariable.hpp
+++ b/SU2_CFD/include/variables/CTurbVariable.hpp
@@ -40,7 +40,7 @@ class CTurbVariable : public CScalarVariable {
   VectorType muT; /*!< \brief Eddy viscosity. */
 
 public:
-  static constexpr size_t MAXNVAR = 2;
+  static constexpr size_t MAXNVAR = 4;
   VectorType turb_index;
   VectorType intermittency;         /*!< \brief Value of the intermittency for the trans. model. */
 
@@ -111,4 +111,19 @@ class CTurbVariable : public CScalarVariable {
    * \param[in] input - Boolean whether In- or Output should be registered.
    */
   void RegisterEddyViscosity(bool input);
+
+    /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   * \param[in] iDim - Dimension index.
+   */
+  inline virtual su2double GetLangevinSourceTermsOld(unsigned long iPoint, unsigned short iDim) const { return 0.0; }
+
+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   * \param[in] iDim - Dimension index.
+   * \param[in] val_stochSource_old - Old value of source term in Langevin equations.
+   */
+  inline virtual void SetLangevinSourceTermsOld(unsigned long iPoint, unsigned short iDim, su2double val_stochSource_old) {}
 };
diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp
index bc4bc55e57c7..a520ed27f465 100644
--- a/SU2_CFD/include/variables/CVariable.hpp
+++ b/SU2_CFD/include/variables/CVariable.hpp
@@ -34,6 +34,7 @@
 #include <cmath>
 #include <iostream>
 #include <cstdlib>
+#include <random>
 
 #include "../../../Common/include/CConfig.hpp"
 #include "../../../Common/include/containers/container_decorators.hpp"
@@ -397,6 +398,46 @@ class CVariable {
    */
   inline virtual void SetDES_LengthScale(unsigned long iPoint, su2double val_des_lengthscale) {}
 
+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   */
+  inline virtual su2double GetLES_Mode(unsigned long iPoint) const { return 0.0; }
+
+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_les_mode - Value of the LES sensor.
+   */
+  inline virtual void SetLES_Mode(unsigned long iPoint, su2double val_les_mode) {}
+
+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   * \param[in] iDim - Dimension index.
+   */
+  inline virtual su2double GetLangevinSourceTerms(unsigned long iPoint, unsigned short iDim) const { return 0.0; }
+
+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   * \param[in] iDim - Dimension index.
+   * \param[in] val_stochSource - Source term in Langevin equations.
+   */
+  inline virtual void SetLangevinSourceTerms(unsigned long iPoint, unsigned short iDim, su2double val_stochSource) {}
+
+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_integral - Value of the integral.
+   */
+  inline virtual void SetBesselIntegral(unsigned long iPoint, su2double val_integral) {}
+
+  /*!
+   * \brief A virtual member.
+   */
+  inline virtual su2double GetBesselIntegral(unsigned long iPoint) const { return 0.0; }
+
   /*!
    * \brief A virtual member.
    * \param[in] iPoint - Point index.
diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp
index 961ce04832ff..f2e1e5b953e8 100644
--- a/SU2_CFD/src/drivers/CDriver.cpp
+++ b/SU2_CFD/src/drivers/CDriver.cpp
@@ -1658,6 +1658,7 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver
           /*--- Incompressible flow, use preconditioning method ---*/
           switch (config->GetKind_Centered_Flow()) {
             case CENTERED::LAX : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentLaxInc_Flow(nDim, nVar_Flow, config); break;
+            case CENTERED::LD2 :
             case CENTERED::JST : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentJSTInc_Flow(nDim, nVar_Flow, config); break;
             default:
               SU2_MPI::Error("Invalid centered scheme or not implemented.\n Currently, only JST and LAX-FRIEDRICH are available for incompressible flows.", CURRENT_FUNCTION);
diff --git a/SU2_CFD/src/numerics/flow/convection/centered.cpp b/SU2_CFD/src/numerics/flow/convection/centered.cpp
index b1ec12d51119..59d904c56050 100644
--- a/SU2_CFD/src/numerics/flow/convection/centered.cpp
+++ b/SU2_CFD/src/numerics/flow/convection/centered.cpp
@@ -311,6 +311,8 @@ CNumerics::ResidualType<> CCentJSTInc_Flow::ComputeResidual(const CConfig* confi
 
   su2double U_i[5] = {0.0}, U_j[5] = {0.0};
   su2double ProjGridVel = 0.0;
+  bool LD2_Scheme = (config->GetKind_Centered_Flow() == CENTERED::LD2);
+  const su2double alpha_LD2 = 0.36;
 
   /*--- Primitive variables at point i and j ---*/
 
@@ -326,6 +328,9 @@ CNumerics::ResidualType<> CCentJSTInc_Flow::ComputeResidual(const CConfig* confi
   for (iDim = 0; iDim < nDim; iDim++) {
     Velocity_i[iDim]    = V_i[iDim+1];
     Velocity_j[iDim]    = V_j[iDim+1];
+  }  
+
+  for (iDim = 0; iDim < nDim; iDim++) {
     MeanVelocity[iDim]  =  0.5*(Velocity_i[iDim]+Velocity_j[iDim]);
     sq_vel_i           += 0.5*Velocity_i[iDim]*Velocity_i[iDim];
     sq_vel_j           += 0.5*Velocity_j[iDim]*Velocity_j[iDim];
@@ -344,6 +349,32 @@ CNumerics::ResidualType<> CCentJSTInc_Flow::ComputeResidual(const CConfig* confi
   MeanCp          = 0.5*(Cp_i          + Cp_j);
   MeanTemperature = 0.5*(Temperature_i + Temperature_j);
 
+  if (LD2_Scheme) {
+    su2double d_ij[3] = {0.0};
+    for (iDim = 0; iDim < nDim; iDim++)
+      d_ij[iDim] = Coord_j[iDim]-Coord_i[iDim];
+    su2double velGrad_i[3][3] = {{0.0}};
+    su2double velGrad_j[3][3] = {{0.0}};
+    su2double pressGrad_i[3] = {0.0};
+    su2double pressGrad_j[3] = {0.0};
+    for (unsigned short jDim = 0; jDim < nDim; jDim++) {
+      pressGrad_i[jDim] = PrimVar_Grad_i[0][jDim];
+      pressGrad_j[jDim] = PrimVar_Grad_j[0][jDim];
+      for (iDim = 0; iDim < nDim; iDim++) {
+        velGrad_i[iDim][jDim] = PrimVar_Grad_i[iDim+1][jDim];
+        velGrad_j[iDim][jDim] = PrimVar_Grad_j[iDim+1][jDim];
+      }
+    }
+    for (iDim = 0; iDim < nDim; iDim++) {
+      MeanVelocity[iDim] += 0.5 * alpha_LD2 * ((velGrad_i[iDim][0] - velGrad_j[iDim][0]) * d_ij[0] +
+                                               (velGrad_i[iDim][1] - velGrad_j[iDim][1]) * d_ij[1] +
+                                               (velGrad_i[iDim][2] - velGrad_j[iDim][2]) * d_ij[2]);
+    }
+    MeanPressure += 0.5 * alpha_LD2 * ((pressGrad_i[0] - pressGrad_j[0]) * d_ij[0] +
+                                       (pressGrad_i[1] - pressGrad_j[1]) * d_ij[1] +
+                                       (pressGrad_i[2] - pressGrad_j[2]) * d_ij[2]);
+  }
+
   /*--- We need the derivative of the equation of state to build the
    preconditioning matrix. For now, the only option is the ideal gas
    law, but in the future, dRhodT should be in the fluid model. ---*/
diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp
index d6bc783ab928..15ba9809b375 100644
--- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp
+++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp
@@ -120,7 +120,8 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar,
                            const su2double* const *val_gradprimvar,
                            const su2double val_turb_ke,
                            const su2double val_laminar_viscosity,
-                           const su2double val_eddy_viscosity) {
+                           const su2double val_eddy_viscosity,
+                           const CConfig* config) {
 
   const su2double Density = val_primvar[nDim+2];
 
@@ -140,6 +141,16 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar,
     // turb_ke is not considered in the stress tensor, see #797
     ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, su2double(0.0));
   }
+
+  /* --- If the Stochastic Backscatter Model is active, add random contribution to stress tensor ---*/
+
+  if (config->GetStochastic_Backscatter()) {
+    for (unsigned short iDim = 0 ; iDim < nDim; iDim++)
+      for (unsigned short jDim = 0 ; jDim < nDim; jDim++) {
+        tau[iDim][jDim] += stochReynStress[iDim][jDim];
+      }
+  }
+
 }
 
 void CAvgGrad_Base::SetHeatFluxVector(const su2double* const* val_gradprimvar, const su2double val_eddy_viscosity,
@@ -445,10 +456,23 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config)
                         Mean_turb_ke, MeanPerturbedRSM);
   }
 
+  /* --- If the Stochastic Backscatter Model is active, add random contribution to stress tensor ---*/
+
+
+  unsigned long timeIter = config->GetTimeIter();
+  unsigned long restartIter = config->GetRestart_Iter();
+  if (config->GetStochastic_Backscatter()) {
+    for (iDim = 0; iDim < nDim; iDim++)
+      Mean_StochVar[iDim] = 0.5*(stochVar_i[iDim] + stochVar_j[iDim]);
+    su2double intensityCoeff = ComputeStochRelaxFactor(config);
+    ComputeStochReynStress(nDim, Mean_PrimVar[nDim+2], Mean_turb_ke,
+                           Mean_StochVar, intensityCoeff, stochReynStress);
+  }
+
   /*--- Get projected flux tensor (viscous residual) ---*/
 
   SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke,
-                  Mean_Laminar_Viscosity, Mean_Eddy_Viscosity);
+                  Mean_Laminar_Viscosity, Mean_Eddy_Viscosity,config);
   if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau);
   if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall);
 
@@ -619,9 +643,23 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi
                         Mean_turb_ke, MeanPerturbedRSM);
   }
 
+  /* --- If the Stochastic Backscatter Model is active, add random contribution to stress tensor ---*/
+
+
+  unsigned long timeIter = config->GetTimeIter();
+  unsigned long restartIter = config->GetRestart_Iter();
+  if (config->GetStochastic_Backscatter()) {
+    for (iDim = 0; iDim < nDim; iDim++)
+      Mean_StochVar[iDim] = 0.5*(stochVar_i[iDim] + stochVar_j[iDim]);
+    su2double intensityCoeff = ComputeStochRelaxFactor(config);
+    ComputeStochReynStress(nDim, Mean_PrimVar[nDim+2], Mean_turb_ke,
+                           Mean_StochVar, intensityCoeff, stochReynStress);
+  }
+
   /*--- Get projected flux tensor (viscous residual) ---*/
+
   SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke,
-                  Mean_Laminar_Viscosity, Mean_Eddy_Viscosity);
+                  Mean_Laminar_Viscosity, Mean_Eddy_Viscosity,config);
   if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau);
   if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall);
 
@@ -944,10 +982,22 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c
                         Mean_turb_ke, MeanPerturbedRSM);
   }
 
+  /* --- If the Stochastic Backscatter Model is active, add random contribution to stress tensor ---*/
+
+  unsigned long timeIter = config->GetTimeIter();
+  unsigned long restartIter = config->GetRestart_Iter();
+  if (config->GetStochastic_Backscatter()) {
+    for (iDim = 0; iDim < nDim; iDim++)
+      Mean_StochVar[iDim] = 0.5*(stochVar_i[iDim] + stochVar_j[iDim]);
+    su2double intensityCoeff = ComputeStochRelaxFactor(config);
+    ComputeStochReynStress(nDim, Mean_PrimVar[nDim+2], Mean_turb_ke,
+                           Mean_StochVar, intensityCoeff, stochReynStress);
+  }
+
   /*--- Get projected flux tensor (viscous residual) ---*/
 
   SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke,
-                  Mean_Laminar_Viscosity, Mean_Eddy_Viscosity);
+                  Mean_Laminar_Viscosity, Mean_Eddy_Viscosity,config);
   if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau);
   if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall);
 
diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp
index d9216b7bdc2f..81573093d745 100644
--- a/SU2_CFD/src/output/CFlowOutput.cpp
+++ b/SU2_CFD/src/output/CFlowOutput.cpp
@@ -1032,6 +1032,14 @@ void CFlowOutput::AddHistoryOutputFields_ScalarRMS_RES(const CConfig* config) {
     case TURB_FAMILY::SA:
       /// DESCRIPTION: Root-mean square residual of nu tilde (SA model).
       AddHistoryOutput("RMS_NU_TILDE", "rms[nu]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of nu tilde (SA model).", HistoryFieldType::RESIDUAL);
+      if (config->GetStochastic_Backscatter()) {
+        /// DESCRIPTION: Root-mean square residual of stochastic vector x-component (Stochastic Backscatter Model).
+        AddHistoryOutput("RMS_STOCH_VAR-X", "rms[stoch_x]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of stochastic vector x-component (Stochastic Backscatter Model).", HistoryFieldType::RESIDUAL);
+        /// DESCRIPTION: Root-mean square residual of stochastic vector y-component (Stochastic Backscatter Model).
+        AddHistoryOutput("RMS_STOCH_VAR-Y", "rms[stoch_y]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of stochastic vector y-component (Stochastic Backscatter Model).", HistoryFieldType::RESIDUAL);
+        /// DESCRIPTION: Root-mean square residual of stochastic vector z-component (Stochastic Backscatter Model).
+        if (nDim==3) AddHistoryOutput("RMS_STOCH_VAR-Z", "rms[stoch_z]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of stochastic vector z-component (Stochastic Backscatter Model).", HistoryFieldType::RESIDUAL);
+      }
       break;
 
     case TURB_FAMILY::KW:
@@ -1087,6 +1095,14 @@ void CFlowOutput::AddHistoryOutputFields_ScalarMAX_RES(const CConfig* config) {
     case TURB_FAMILY::SA:
       /// DESCRIPTION: Maximum residual of nu tilde (SA model).
       AddHistoryOutput("MAX_NU_TILDE", "max[nu]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of nu tilde (SA model).", HistoryFieldType::RESIDUAL);
+      if (config->GetStochastic_Backscatter()) {
+        /// DESCRIPTION: Maximum residual of stochastic vector x-component (Stochastic Backscatter Model).
+        AddHistoryOutput("MAX_STOCH_VAR-X", "max[stoch_x]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of stochastic vector x-component (Stochastic Backscatter Model).", HistoryFieldType::RESIDUAL);
+        /// DESCRIPTION: Maximum residual of stochastic vector y-component (Stochastic Backscatter Model).
+        AddHistoryOutput("MAX_STOCH_VAR-Y", "max[stoch_y]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of stochastic vector y-component (Stochastic Backscatter Model).", HistoryFieldType::RESIDUAL);
+        /// DESCRIPTION: Maximum residual of stochastic vector z-component (Stochastic Backscatter Model).
+        if (nDim==3) AddHistoryOutput("MAX_STOCH_VAR-Z", "max[stoch_z]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of stochastic vector z-component (Stochastic Backscatter Model).", HistoryFieldType::RESIDUAL);
+      }
       break;
 
     case TURB_FAMILY::KW:
@@ -1228,8 +1244,21 @@ void CFlowOutput::LoadHistoryDataScalar(const CConfig* config, const CSolver* co
     case TURB_FAMILY::SA:
       SetHistoryOutputValue("RMS_NU_TILDE", log10(solver[TURB_SOL]->GetRes_RMS(0)));
       SetHistoryOutputValue("MAX_NU_TILDE", log10(solver[TURB_SOL]->GetRes_Max(0)));
+      if (config->GetStochastic_Backscatter()) {
+        SetHistoryOutputValue("RMS_STOCH_VAR-X", log10(solver[TURB_SOL]->GetRes_RMS(1)));
+        SetHistoryOutputValue("RMS_STOCH_VAR-Y", log10(solver[TURB_SOL]->GetRes_RMS(2)));
+        if (nDim==3) SetHistoryOutputValue("RMS_STOCH_VAR-Z", log10(solver[TURB_SOL]->GetRes_RMS(3)));
+        SetHistoryOutputValue("MAX_STOCH_VAR-X", log10(solver[TURB_SOL]->GetRes_Max(1)));
+        SetHistoryOutputValue("MAX_STOCH_VAR-Y", log10(solver[TURB_SOL]->GetRes_Max(2)));
+        if (nDim==3) SetHistoryOutputValue("MAX_STOCH_VAR-Z", log10(solver[TURB_SOL]->GetRes_Max(3)));
+      }
       if (multiZone) {
         SetHistoryOutputValue("BGS_NU_TILDE", log10(solver[TURB_SOL]->GetRes_BGS(0)));
+        if (config->GetStochastic_Backscatter()) {
+          SetHistoryOutputValue("BGS_STOCH_VAR-X", log10(solver[TURB_SOL]->GetRes_BGS(1)));
+          SetHistoryOutputValue("BGS_STOCH_VAR-Y", log10(solver[TURB_SOL]->GetRes_BGS(2)));
+          if (nDim==3) SetHistoryOutputValue("BGS_STOCH_VAR-Z", log10(solver[TURB_SOL]->GetRes_BGS(3)));
+        }
       }
       break;
 
@@ -1551,6 +1580,16 @@ void CFlowOutput::SetVolumeOutputFieldsScalarMisc(const CConfig* config) {
   if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES) {
     AddVolumeOutput("DES_LENGTHSCALE", "DES_LengthScale", "DDES", "DES length scale value");
     AddVolumeOutput("WALL_DISTANCE", "Wall_Distance", "DDES", "Wall distance value");
+    AddVolumeOutput("LES_SENSOR","LES_Sensor","DDES","LES sensor value");
+    if (config->GetStochastic_Backscatter()) {
+      AddVolumeOutput("STOCHVAR_X", "StochVar_x", "BACKSCATTER", "x-component of the stochastic vector potential");
+      AddVolumeOutput("STOCHVAR_Y", "StochVar_y", "BACKSCATTER", "y-component of the stochastic vector potential");
+      if (nDim==3) AddVolumeOutput("STOCHVAR_Z", "StochVar_z", "BACKSCATTER", "z-component of the stochastic vector potential");
+      AddVolumeOutput("STOCHSOURCE_X", "StochSource_x", "BACKSCATTER", "x-component of the stochastic source vector");
+      AddVolumeOutput("STOCHSOURCE_Y", "StochSource_y", "BACKSCATTER", "y-component of the stochastic source vector");
+      if (nDim==3) AddVolumeOutput("STOCHSOURCE_Z", "StochSource_z", "BACKSCATTER", "z-component of the stochastic source vector");
+      AddVolumeOutput("ENERGY_BACKSCATTER_RATIO", "Energy_Backscatter_Ratio", "BACKSCATTER", "Energy backscatter from unresolved to resolved scales (divided by the turbulent dissipation of resolved kinetic energy)");
+    }
   }
 
   if (config->GetViscous()) {
@@ -1650,6 +1689,93 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
   if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES) {
     SetVolumeOutputValue("DES_LENGTHSCALE", iPoint, Node_Flow->GetDES_LengthScale(iPoint));
     SetVolumeOutputValue("WALL_DISTANCE", iPoint, Node_Geo->GetWall_Distance(iPoint));
+    SetVolumeOutputValue("LES_SENSOR", iPoint, Node_Flow->GetLES_Mode(iPoint));
+    if (config->GetStochastic_Backscatter()) {
+      SetVolumeOutputValue("STOCHVAR_X", iPoint, Node_Turb->GetSolution(iPoint, 1));
+      SetVolumeOutputValue("STOCHVAR_Y", iPoint, Node_Turb->GetSolution(iPoint, 2));
+      if (nDim==3) SetVolumeOutputValue("STOCHVAR_Z", iPoint, Node_Turb->GetSolution(iPoint, 3));
+      SetVolumeOutputValue("STOCHSOURCE_X", iPoint, Node_Turb->GetLangevinSourceTerms(iPoint, 0));
+      SetVolumeOutputValue("STOCHSOURCE_Y", iPoint, Node_Turb->GetLangevinSourceTerms(iPoint, 1));
+      if (nDim==3) SetVolumeOutputValue("STOCHSOURCE_Z", iPoint, Node_Turb->GetLangevinSourceTerms(iPoint, 2));
+      const su2double rho = Node_Flow->GetDensity(iPoint);
+      const su2double nu_t = Node_Flow->GetEddyViscosity(iPoint) / rho;
+      const su2double DES_lengthscale = max(Node_Flow->GetDES_LengthScale(iPoint), 1e-10);
+      const su2double lesSensor = Node_Flow->GetLES_Mode(iPoint);
+      const su2double mag = config->GetSBS_Cmag();
+      const su2double threshold = config->GetStochFdThreshold();
+      const auto VelocityGradient = Node_Flow->GetVelocityGradient(iPoint);
+      su2double strainMag2 = 0.0;
+      for (unsigned long iDim = 0; iDim < nDim; iDim++) {
+        strainMag2 += pow(VelocityGradient(iDim, iDim), 2);
+      }
+      strainMag2 += 2.0*pow(0.5*(VelocityGradient(0,1) + VelocityGradient(1,0)), 2);
+      if (nDim == 3) {
+        strainMag2 += 2.0*pow(0.5*(VelocityGradient(0,2) + VelocityGradient(2,0)), 2);
+        strainMag2 += 2.0*pow(0.5*(VelocityGradient(1,2) + VelocityGradient(2,1)), 2);
+      }
+      strainMag2 *= 2.0;
+      su2double tke_estim = 0.0;
+      if (lesSensor > threshold) tke_estim = sqrt(strainMag2) * nu_t;
+      const su2double csi_x = Node_Turb->GetSolution(iPoint, 1);
+      const su2double csi_y = Node_Turb->GetSolution(iPoint, 2);
+      const su2double csi_z = Node_Turb->GetSolution(iPoint, 3);
+      const su2double R_xy =   mag * tke_estim * csi_z;
+      const su2double R_xz = - mag * tke_estim * csi_y;
+      const su2double R_yz =   mag * tke_estim * csi_x;
+      const su2double energy_res_to_mod = nu_t * strainMag2;
+      const auto vorticity = Node_Flow->GetVorticity(iPoint);
+      const su2double energy_backscatter = R_xy*vorticity[2] - R_xz*vorticity[1] + R_yz*vorticity[0];
+      const su2double energy_backscatter_ratio = energy_backscatter / (energy_res_to_mod + 1e-10);
+      SetVolumeOutputValue("ENERGY_BACKSCATTER_RATIO", iPoint, energy_backscatter_ratio);
+    }
+  }
+
+  if (config->GetTime_Domain()) {  
+    const su2double rho = Node_Flow->GetDensity(iPoint);
+    const su2double nu_t = Node_Flow->GetEddyViscosity(iPoint) / rho;
+    const auto vel_grad = Node_Flow->GetVelocityGradient(iPoint);
+    const su2double vel_div = vel_grad(0,0) + vel_grad(1,1) + (nDim ==3 ? vel_grad(2,2) : 0.0);
+    const su2double tau_xx = nu_t * (2*vel_grad(0,0) - (2.0/3.0)*vel_div);
+    const su2double tau_yy = nu_t * (2*vel_grad(1,1) - (2.0/3.0)*vel_div);
+    const su2double tau_xy = nu_t * (vel_grad(0,1) + vel_grad(1,0));
+    SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_11", iPoint, -tau_xx);
+    SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_22", iPoint, -tau_yy);
+    SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_12", iPoint, -tau_xy);
+    if (nDim == 3){
+      const su2double tau_zz = nu_t * (2*vel_grad(2,2) - (2.0/3.0)*vel_div);
+      const su2double tau_xz = nu_t * (vel_grad(0,2) + vel_grad(2,0));
+      const su2double tau_yz = nu_t * (vel_grad(1,2) + vel_grad(2,1));
+      SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_33", iPoint, -tau_zz);
+      SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_13", iPoint, -tau_xz);
+      SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_23", iPoint, -tau_yz);
+    }
+    if (config->GetKind_HybridRANSLES()!=NO_HYBRIDRANSLES && config->GetStochastic_Backscatter()) {
+      const su2double DES_lengthscale = max(Node_Flow->GetDES_LengthScale(iPoint), 1e-10);
+      const su2double lesSensor = Node_Flow->GetLES_Mode(iPoint);
+      const su2double mag = config->GetSBS_Cmag();
+      const su2double threshold = config->GetStochFdThreshold();
+      su2double tke_estim = 0.0;
+      su2double strainMag2 = 0.0;
+      for (unsigned long iDim = 0; iDim < nDim; iDim++) {
+        strainMag2 += pow(vel_grad(iDim, iDim), 2);
+      }
+      strainMag2 += 2.0*pow(0.5*(vel_grad(0,1) + vel_grad(1,0)), 2);
+      if (nDim == 3) {
+        strainMag2 += 2.0*pow(0.5*(vel_grad(0,2) + vel_grad(2,0)), 2);
+        strainMag2 += 2.0*pow(0.5*(vel_grad(1,2) + vel_grad(2,1)), 2);
+      }
+      strainMag2 *= 2.0;
+      if (lesSensor > threshold) tke_estim = sqrt(strainMag2) * nu_t;
+      const su2double csi_x = Node_Turb->GetSolution(iPoint, 1);
+      const su2double csi_y = Node_Turb->GetSolution(iPoint, 2);
+      const su2double csi_z = Node_Turb->GetSolution(iPoint, 3);
+      const su2double R_xy =   mag * tke_estim * csi_z;
+      const su2double R_xz = - mag * tke_estim * csi_y;
+      const su2double R_yz =   mag * tke_estim * csi_x;
+      SetAvgVolumeOutputValue("STOCHASTIC_REYNOLDS_STRESS_12", iPoint, -R_xy);
+      SetAvgVolumeOutputValue("STOCHASTIC_REYNOLDS_STRESS_13", iPoint, -R_xz);
+      SetAvgVolumeOutputValue("STOCHASTIC_REYNOLDS_STRESS_23", iPoint, -R_yz);
+    }
   }
 
   switch (config->GetKind_Species_Model()) {
@@ -1729,6 +1855,13 @@ void CFlowOutput::LoadSurfaceData(CConfig *config, CGeometry *geometry, CSolver
     SetVolumeOutputValue("SKIN_FRICTION-Z", iPoint, solver[FLOW_SOL]->GetCSkinFriction(iMarker, iVertex, 2));
   SetVolumeOutputValue("HEAT_FLUX", iPoint, solver[heat_sol]->GetHeatFlux(iMarker, iVertex));
   SetVolumeOutputValue("Y_PLUS", iPoint, solver[FLOW_SOL]->GetYPlus(iMarker, iVertex));
+  
+  if (config->GetTime_Domain()) {
+    SetAvgVolumeOutputValue("MEAN_SKIN_FRICTION-X", iPoint, solver[FLOW_SOL]->GetCSkinFriction(iMarker, iVertex, 0));
+    SetAvgVolumeOutputValue("MEAN_SKIN_FRICTION-Y", iPoint, solver[FLOW_SOL]->GetCSkinFriction(iMarker, iVertex, 1));
+    if (nDim == 3) 
+      SetAvgVolumeOutputValue("MEAN_SKIN_FRICTION-Z", iPoint, solver[FLOW_SOL]->GetCSkinFriction(iMarker, iVertex, 2));
+  }
 }
 
 void CFlowOutput::AddAerodynamicCoefficients(const CConfig* config) {
@@ -4036,6 +4169,11 @@ void CFlowOutput::SetTimeAveragedFields() {
     AddVolumeOutput("MEAN_VELOCITY-Z", "MeanVelocity_z", "TIME_AVERAGE", "Mean velocity z-component");
 
   AddVolumeOutput("MEAN_PRESSURE", "MeanPressure", "TIME_AVERAGE", "Mean pressure");
+  AddVolumeOutput("MEAN_SKIN_FRICTION-X", "MeanSkinFriction_x", "TIME_AVERAGE", "Mean skin friction x-component");
+  AddVolumeOutput("MEAN_SKIN_FRICTION-Y", "MeanSkinFriction_y", "TIME_AVERAGE", "Mean skin friction y-component");
+  if (nDim==3)
+    AddVolumeOutput("MEAN_SKIN_FRICTION-Z", "MeanSkinFriction_z", "TIME_AVERAGE", "Mean skin friction z-component");
+
   AddVolumeOutput("RMS_U",   "RMS[u]", "TIME_AVERAGE", "RMS u");
   AddVolumeOutput("RMS_V",   "RMS[v]", "TIME_AVERAGE", "RMS v");
   AddVolumeOutput("RMS_UV",  "RMS[uv]", "TIME_AVERAGE", "RMS uv");
@@ -4052,6 +4190,19 @@ void CFlowOutput::SetTimeAveragedFields() {
     AddVolumeOutput("UWPRIME", "w'u'", "TIME_AVERAGE", "Mean Reynolds-stress component w'u'");
     AddVolumeOutput("VWPRIME", "w'v'", "TIME_AVERAGE", "Mean Reynolds-stress component w'v'");
   }
+
+  AddVolumeOutput("MODELED_REYNOLDS_STRESS_11", "ModeledReynoldsStress_11", "TIME_AVERAGE", "Modeled Reynolds stress xx-component");
+  AddVolumeOutput("MODELED_REYNOLDS_STRESS_22", "ModeledReynoldsStress_22", "TIME_AVERAGE", "Modeled Reynolds stress yy-component");
+  AddVolumeOutput("MODELED_REYNOLDS_STRESS_12", "ModeledReynoldsStress_12", "TIME_AVERAGE", "Modeled Reynolds stress xy-component");
+  if (nDim == 3){
+    AddVolumeOutput("MODELED_REYNOLDS_STRESS_33", "ModeledReynoldsStress_33", "TIME_AVERAGE", "Modeled Reynolds stress zz-component");
+    AddVolumeOutput("MODELED_REYNOLDS_STRESS_13", "ModeledReynoldsStress_13", "TIME_AVERAGE", "Modeled Reynolds stress xz-component");
+    AddVolumeOutput("MODELED_REYNOLDS_STRESS_23", "ModeledReynoldsStress_23", "TIME_AVERAGE", "Modeled Reynolds stress yz-component");
+  }
+
+  AddVolumeOutput("STOCHASTIC_REYNOLDS_STRESS_12", "StochasticReynoldsStress_12", "TIME_AVERAGE", "Stochastic Reynolds stress xy-component");
+  AddVolumeOutput("STOCHASTIC_REYNOLDS_STRESS_13", "StochasticReynoldsStress_13", "TIME_AVERAGE", "Stochastic Reynolds stress xz-component");
+  AddVolumeOutput("STOCHASTIC_REYNOLDS_STRESS_23", "StochasticReynoldsStress_23", "TIME_AVERAGE", "Stochastic Reynolds stress yz-component");
 }
 
 void CFlowOutput::LoadTimeAveragedData(unsigned long iPoint, const CVariable *Node_Flow){
@@ -4062,7 +4213,6 @@ void CFlowOutput::LoadTimeAveragedData(unsigned long iPoint, const CVariable *No
     SetAvgVolumeOutputValue("MEAN_VELOCITY-Z", iPoint, Node_Flow->GetVelocity(iPoint,2));
 
   SetAvgVolumeOutputValue("MEAN_PRESSURE", iPoint, Node_Flow->GetPressure(iPoint));
-
   SetAvgVolumeOutputValue("RMS_U", iPoint, pow(Node_Flow->GetVelocity(iPoint,0),2));
   SetAvgVolumeOutputValue("RMS_V", iPoint, pow(Node_Flow->GetVelocity(iPoint,1),2));
   SetAvgVolumeOutputValue("RMS_UV", iPoint, Node_Flow->GetVelocity(iPoint,0) * Node_Flow->GetVelocity(iPoint,1));
diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
index c99f87896c2d..f771594113de 100644
--- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
@@ -959,7 +959,7 @@ void CIncEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_
 
   const bool implicit   = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
   const bool center     = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED);
-  const bool center_jst = (config->GetKind_Centered_Flow() == CENTERED::JST) && (iMesh == MESH_0);
+  const bool center_jst = (config->GetKind_Centered_Flow() == CENTERED::JST || config->GetKind_Centered_Flow() == CENTERED::LD2) && (iMesh == MESH_0);
   const bool outlet     = (config->GetnMarker_Outlet() != 0);
 
   /*--- Set the primitive variables ---*/
@@ -1132,8 +1132,9 @@ void CIncEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_co
   unsigned long iPoint, jPoint;
 
   const bool implicit    = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
-  const bool jst_scheme  = ((config->GetKind_Centered_Flow() == CENTERED::JST) && (iMesh == MESH_0));
+  const bool jst_scheme  = ((config->GetKind_Centered_Flow() == CENTERED::JST || config->GetKind_Centered_Flow() == CENTERED::LD2) && (iMesh == MESH_0));
   const bool bounded_scalar = config->GetBounded_Scalar();
+  const bool LD2_Scheme = (config->GetKind_Centered_Flow() == CENTERED::LD2);
 
   /*--- For hybrid parallel AD, pause preaccumulation if there is shared reading of
   * variables, otherwise switch to the faster adjoint evaluation mode. ---*/
@@ -1171,6 +1172,16 @@ void CIncEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_co
       numerics->SetSensor(nodes->GetSensor(iPoint), nodes->GetSensor(jPoint));
     }
 
+    if (LD2_Scheme) {
+      numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(jPoint));
+      if (!geometry->nodes->GetPeriodicBoundary(iPoint) || (geometry->nodes->GetPeriodicBoundary(iPoint) 
+          && !geometry->nodes->GetPeriodicBoundary(jPoint))) {
+        numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(jPoint));
+      } else {
+        numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint));
+      }
+    }
+
     /*--- Grid movement ---*/
 
     if (dynamic_grid) {
diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp
index 6477fa6f3619..8c2e40f90dcb 100644
--- a/SU2_CFD/src/solvers/CIncNSSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp
@@ -360,7 +360,7 @@ void CIncNSSolver::Compute_Enthalpy_Diffusion(unsigned long iEdge, CGeometry* ge
 unsigned long CIncNSSolver::SetPrimitive_Variables(CSolver **solver_container, const CConfig *config) {
 
   unsigned long iPoint, nonPhysicalPoints = 0;
-  su2double eddy_visc = 0.0, turb_ke = 0.0, DES_LengthScale = 0.0;
+  su2double eddy_visc = 0.0, turb_ke = 0.0, DES_LengthScale = 0.0, LES_Mode = 0.0;
   const su2double* scalar = nullptr;
   const TURB_MODEL turb_model = config->GetKind_Turb_Model();
   const SPECIES_MODEL species_model = config->GetKind_Species_Model();
@@ -380,6 +380,7 @@ unsigned long CIncNSSolver::SetPrimitive_Variables(CSolver **solver_container, c
 
       if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES){
         DES_LengthScale = solver_container[TURB_SOL]->GetNodes()->GetDES_LengthScale(iPoint);
+        LES_Mode = solver_container[TURB_SOL]->GetNodes()->GetLES_Mode(iPoint); 
       }
     }
 
@@ -400,6 +401,10 @@ unsigned long CIncNSSolver::SetPrimitive_Variables(CSolver **solver_container, c
 
     nodes->SetDES_LengthScale(iPoint,DES_LengthScale);
 
+    /*--- Set the LES sensor ---*/
+
+    nodes->SetLES_Mode(iPoint, LES_Mode);
+
   }
   END_SU2_OMP_FOR
 
diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp
index 7a1309027a65..c4ccfc900202 100644
--- a/SU2_CFD/src/solvers/CNSSolver.cpp
+++ b/SU2_CFD/src/solvers/CNSSolver.cpp
@@ -150,6 +150,8 @@ unsigned long CNSSolver::SetPrimitive_Variables(CSolver **solver_container, cons
       if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES) {
         su2double DES_LengthScale = solver_container[TURB_SOL]->GetNodes()->GetDES_LengthScale(iPoint);
         nodes->SetDES_LengthScale(iPoint, DES_LengthScale);
+        su2double LES_Mode = solver_container[TURB_SOL]->GetNodes()->GetLES_Mode(iPoint);
+        nodes->SetLES_Mode(iPoint, LES_Mode);
       }
     }
 
diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp
index f544b694ba82..d4a23d6de3f3 100644
--- a/SU2_CFD/src/solvers/CSolver.cpp
+++ b/SU2_CFD/src/solvers/CSolver.cpp
@@ -363,8 +363,8 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
 
   auto *Diff      = new su2double[nVar];
   auto *Und_Lapl  = new su2double[nVar];
-  auto *Sol_Min   = new su2double[nPrimVarGrad];
-  auto *Sol_Max   = new su2double[nPrimVarGrad];
+  auto *Sol_Min   = new su2double[nVar];
+  auto *Sol_Max   = new su2double[nVar];
   auto *rotPrim_i = new su2double[nPrimVar];
   auto *rotPrim_j = new su2double[nPrimVar];
 
@@ -1356,6 +1356,10 @@ void CSolver::GetCommCountAndType(const CConfig* config,
       COUNT_PER_POINT  = nVar+1;
       MPI_TYPE         = COMM_TYPE_DOUBLE;
       break;
+    case MPI_QUANTITIES::STOCH_SOURCE_LANG:
+      COUNT_PER_POINT  = nDim;
+      MPI_TYPE         = COMM_TYPE_DOUBLE;
+      break;
     case MPI_QUANTITIES::SOLUTION_FEA:
       if (config->GetTime_Domain())
         COUNT_PER_POINT  = nVar*3;
@@ -1483,6 +1487,10 @@ void CSolver::InitiateComms(CGeometry *geometry,
               bufDSend[buf_offset+iVar] = base_nodes->GetSolution(iPoint, iVar);
             bufDSend[buf_offset+nVar]   = base_nodes->GetmuT(iPoint);
             break;
+          case MPI_QUANTITIES::STOCH_SOURCE_LANG:
+            for (iDim = 0; iDim < nDim; iDim++)
+              bufDSend[buf_offset+iDim] = base_nodes->GetLangevinSourceTerms(iPoint, iDim);
+            break;
           case MPI_QUANTITIES::UNDIVIDED_LAPLACIAN:
             for (iVar = 0; iVar < nVar; iVar++)
               bufDSend[buf_offset+iVar] = base_nodes->GetUndivided_Laplacian(iPoint, iVar);
@@ -1631,6 +1639,10 @@ void CSolver::CompleteComms(CGeometry *geometry,
               base_nodes->SetSolution(iPoint, iVar, bufDRecv[buf_offset+iVar]);
             base_nodes->SetmuT(iPoint,bufDRecv[buf_offset+nVar]);
             break;
+          case MPI_QUANTITIES::STOCH_SOURCE_LANG:
+            for (iDim = 0; iDim < nDim; iDim++)
+              base_nodes->SetLangevinSourceTerms(iPoint, iDim, bufDRecv[buf_offset+iDim]);
+            break;
           case MPI_QUANTITIES::UNDIVIDED_LAPLACIAN:
             for (iVar = 0; iVar < nVar; iVar++)
               base_nodes->SetUnd_Lapl(iPoint, iVar, bufDRecv[buf_offset+iVar]);
diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp
index 5ea64882b0b7..b197f8f8a0e1 100644
--- a/SU2_CFD/src/solvers/CTurbSASolver.cpp
+++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp
@@ -30,6 +30,7 @@
 #include "../../include/variables/CFlowVariable.hpp"
 #include "../../../Common/include/parallelization/omp_structure.hpp"
 #include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
+#include "../../../Common/include/toolboxes/random_toolbox.hpp"
 
 
 CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned short iMesh, CFluidModel* FluidModel)
@@ -54,6 +55,16 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor
 
   nDim = geometry->GetnDim();
 
+  /*--- Add Langevin equations if the Stochastic Backscatter Model is used ---*/
+
+  if (config->GetStochastic_Backscatter()) { 
+    if (nDim == 3) {
+      nVar += nDim; nVarGrad = nPrimVar = nVar;
+    } else {
+      SU2_MPI::Error("Stochastic Backscatter Model available for 3D flows only.", CURRENT_FUNCTION);
+    }
+  }
+
   /*--- Single grid simulation ---*/
 
   if (iMesh == MESH_0 || config->GetMGCycle() == FULLMG_CYCLE) {
@@ -109,17 +120,22 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor
   }
 
   Solution_Inf[0] = nu_tilde_Inf;
+  if (config->GetStochastic_Backscatter()) {
+    for (unsigned short iVar = 1; iVar < nVar; iVar++) {
+      Solution_Inf[iVar] = 0.0;
+    }
+  }
 
   /*--- Factor_nu_Engine ---*/
-  Factor_nu_Engine = config->GetNuFactor_Engine();
-  nu_tilde_Engine  = Factor_nu_Engine*Viscosity_Inf/Density_Inf;
+  Factor_nu_Engine   = config->GetNuFactor_Engine();
+  nu_tilde_Engine[0] = Factor_nu_Engine*Viscosity_Inf/Density_Inf;
   if (config->GetSAParsedOptions().bc) {
-    nu_tilde_Engine  = 0.005*Factor_nu_Engine*Viscosity_Inf/Density_Inf;
+    nu_tilde_Engine[0] = 0.005*Factor_nu_Engine*Viscosity_Inf/Density_Inf;
   }
 
   /*--- Factor_nu_ActDisk ---*/
-  Factor_nu_ActDisk = config->GetNuFactor_Engine();
-  nu_tilde_ActDisk  = Factor_nu_ActDisk*Viscosity_Inf/Density_Inf;
+  Factor_nu_ActDisk   = config->GetNuFactor_Engine();
+  nu_tilde_ActDisk[0] = Factor_nu_ActDisk*Viscosity_Inf/Density_Inf;
 
   /*--- Eddy viscosity at infinity ---*/
   su2double Ji, Ji_3, fv1, cv1_3 = 7.1*7.1*7.1;
@@ -155,8 +171,16 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor
    * due to arbitrary number of turbulence variables ---*/
 
   Inlet_TurbVars.resize(nMarker);
-  for (unsigned long iMarker = 0; iMarker < nMarker; iMarker++)
+  for (unsigned long iMarker = 0; iMarker < nMarker; iMarker++) {
     Inlet_TurbVars[iMarker].resize(nVertex[iMarker],nVar) = nu_tilde_Inf;
+    if (config->GetStochastic_Backscatter()) {
+      for (unsigned long iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) {
+        for (unsigned short iVar = 1; iVar < nVar; iVar++) {
+          Inlet_TurbVars[iMarker](iVertex,iVar) = 0.0;
+        }
+      }
+    }
+  }
 
   /*--- Store the initial CFL number for all grid points. ---*/
 
@@ -203,6 +227,34 @@ void CTurbSASolver::Preprocessing(CGeometry *geometry, CSolver **solver_containe
 
     SetDES_LengthScale(solver_container, geometry, config);
 
+    /*--- Compute source terms for Langevin equations ---*/
+
+    bool backscatter = config->GetStochastic_Backscatter();
+    unsigned long innerIter = config->GetInnerIter();
+    unsigned long timeIter = config->GetTimeIter();
+    unsigned long restartIter = config->GetRestart_Iter();
+
+    if (backscatter && innerIter==0) {
+      SetLangevinSourceTerms(config, geometry);
+      const unsigned short maxIter = config->GetSBS_maxIterSmooth();
+      bool dual_time = ((config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) ||
+                        (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND));
+      if (maxIter>0) SmoothLangevinSourceTerms(config, geometry);
+      if (timeIter == restartIter) {
+        for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
+          for (unsigned short iVar = 1; iVar < nVar; iVar++) {
+            const su2double randomSource = nodes->GetLangevinSourceTerms(iPoint, iVar-1);
+            nodes->SetSolution(iPoint, iVar, randomSource);
+            nodes->SetSolution_Old(iPoint, iVar, randomSource);
+            if (dual_time) {
+              nodes->Set_Solution_time_n(iPoint, iVar, randomSource);
+              nodes->Set_Solution_time_n1(iPoint, iVar, randomSource);
+            }
+          }
+        }
+      }
+    }
+
   }
 
 }
@@ -378,6 +430,14 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai
 
       numerics->SetDistance(nodes->GetDES_LengthScale(iPoint), 0.0);
 
+      /*--- Compute source terms in Langevin equations (Stochastic Basckscatter Model) ---*/
+
+      if (config->GetStochastic_Backscatter()) {
+        for (unsigned short iDim = 0; iDim < nDim; iDim++)
+          numerics->SetStochSource(nodes->GetLangevinSourceTerms(iPoint, iDim), iDim);
+        numerics->SetLES_Mode(nodes->GetLES_Mode(iPoint), 0.0);
+      }
+
     }
 
     /*--- Effective Intermittency ---*/
@@ -401,7 +461,7 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai
     if (transition_BC || config->GetKind_Trans_Model() != TURB_TRANS_MODEL::NONE) {
       nodes->SetIntermittency(iPoint,numerics->GetIntermittencyEff());
     }
-
+    
     /*--- Subtract residual and the Jacobian ---*/
 
     LinSysRes.SubtractBlock(iPoint, residual);
@@ -503,19 +563,21 @@ void CTurbSASolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_conta
          su2double coeff = (nu_total/sigma);
          su2double RoughWallBC = nodes->GetSolution(iPoint,0)/(0.03*Roughness_Height);
 
-         su2double Res_Wall;// = new su2double [nVar];
-         Res_Wall = coeff*RoughWallBC*Area;
-         LinSysRes.SubtractBlock(iPoint, &Res_Wall);
+         su2double Res_Wall[MAXNVAR] = {0.0};
+         Res_Wall[0] = coeff*RoughWallBC*Area;
+         LinSysRes.SubtractBlock(iPoint, Res_Wall);
 
-         su2double Jacobian_i = (laminar_viscosity /density *Area)/(0.03*Roughness_Height*sigma);
-         Jacobian_i += 2.0*RoughWallBC*Area/sigma;
-         if (implicit) Jacobian.AddVal2Diag(iPoint, -Jacobian_i);
+         Jacobian_i[0][0] = (laminar_viscosity /density *Area)/(0.03*Roughness_Height*sigma);
+         Jacobian_i[0][0] += 2.0*RoughWallBC*Area/sigma;
+         if (implicit) Jacobian.AddVal2Diag(iPoint, 0, -Jacobian_i[0][0]);
       }
     }
   }
   END_SU2_OMP_FOR
 }
 
+
+
 void CTurbSASolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics,
                                        CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) {
 
@@ -558,7 +620,7 @@ void CTurbSASolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CN
       conv_numerics->SetPrimitive(V_domain, V_inlet);
 
       /*--- Non-dimensionalize Inlet_TurbVars if Inlet-Files are used. ---*/
-      su2double Inlet_Vars[MAXNVAR];
+      su2double Inlet_Vars[MAXNVAR] = {0.0};
       Inlet_Vars[0] = Inlet_TurbVars[val_marker][iVertex][0];
       if (config->GetInlet_Profile_From_File()) {
          Inlet_Vars[0] *= config->GetDensity_Ref() / config->GetViscosity_Ref();
@@ -863,7 +925,7 @@ void CTurbSASolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_cont
 
       /*--- Set the turbulent variable states (prescribed for an inflow) ---*/
 
-      conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), &nu_tilde_Engine);
+      conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), nu_tilde_Engine);
 
       /*--- Set various other quantities in the conv_numerics class ---*/
 
@@ -1016,7 +1078,7 @@ void CTurbSASolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container,
     }
     else {
       /*--- Outflow analysis ---*/
-      conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), &nu_tilde_ActDisk);
+      conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), nu_tilde_ActDisk);
     }
 
     /*--- Grid Movement ---*/
@@ -1075,7 +1137,8 @@ void CTurbSASolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_c
   /*--- Loop over all the vertices on this boundary marker ---*/
   for (auto iSpan = 0u; iSpan < nSpanWiseSections ; iSpan++){
 
-    su2double extAverageNu = solver_container[FLOW_SOL]->GetExtAverageNu(val_marker, iSpan);
+    su2double extAverageNu[MAXNVAR] = {0.0};
+    extAverageNu[0] = solver_container[FLOW_SOL]->GetExtAverageNu(val_marker, iSpan);
 
     /*--- Loop over all the vertices on this boundary marker ---*/
 
@@ -1111,7 +1174,7 @@ void CTurbSASolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_c
 
       /*--- Set the turbulent variable states (prescribed for an inflow) ---*/
 
-      conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), &extAverageNu);
+      conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), extAverageNu);
 
       /*--- Set various other quantities in the conv_numerics class ---*/
 
@@ -1141,7 +1204,7 @@ void CTurbSASolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_c
 
       /*--- Turbulent variables w/o reconstruction, and its gradients ---*/
 
-      visc_numerics->SetScalarVar(nodes->GetSolution(iPoint), &extAverageNu);
+      visc_numerics->SetScalarVar(nodes->GetSolution(iPoint), extAverageNu);
 
       visc_numerics->SetScalarVarGradient(nodes->GetGradient(iPoint),
                                         nodes->GetGradient(iPoint));
@@ -1180,7 +1243,8 @@ void CTurbSASolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_contain
     FluidModel->SetTDState_Prho(pressure, rho);
     su2double muLam = FluidModel->GetLaminarViscosity();
 
-    su2double nu_tilde  = Factor_nu_Inf*muLam/rho;
+    su2double nu_tilde[MAXNVAR] = {0.0};
+    nu_tilde[0] = Factor_nu_Inf*muLam/rho;
 
     /*--- Loop over all the vertices on this boundary marker ---*/
 
@@ -1216,7 +1280,7 @@ void CTurbSASolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_contain
 
       /*--- Set the turbulent variable states (prescribed for an inflow) ---*/
 
-      conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), &nu_tilde);
+      conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), nu_tilde);
 
       if (dynamic_grid)
         conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint),
@@ -1246,7 +1310,7 @@ void CTurbSASolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_contain
 
       /*--- Turbulent variables w/o reconstruction, and its gradients ---*/
 
-      visc_numerics->SetScalarVar(nodes->GetSolution(iPoint), &nu_tilde);
+      visc_numerics->SetScalarVar(nodes->GetSolution(iPoint), nu_tilde);
 
       visc_numerics->SetScalarVarGradient(nodes->GetGradient(iPoint),
                                           nodes->GetGradient(iPoint));
@@ -1337,7 +1401,9 @@ void CTurbSASolver::SetTurbVars_WF(CGeometry *geometry, CSolver **solver_contain
         if (counter > max_iter) break;
       }
 
-      nodes->SetSolution_Old(iPoint_Neighbor, &nu_til);
+      su2double nuTil[MAXNVAR] = {0.0};
+      nuTil[0] = nu_til;
+      nodes->SetSolution_Old(iPoint_Neighbor, nuTil);
       LinSysRes.SetBlock_Zero(iPoint_Neighbor);
 
       /*--- includes 1 in the diagonal ---*/
@@ -1396,7 +1462,9 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
     su2double psi_2 = (1.0 - (cb1/(cw1*k2*fw_star))*(ft2 + (1.0 - ft2)*fv2))/(fv1 * max(1.0e-10,1.0-ft2));
     psi_2 = min(100.0,psi_2);
 
-    su2double lengthScale = 0.0;
+    su2double lengthScale = 0.0, lesSensor = 0.0;
+
+    const su2double LES_FilterWidth = config->GetLES_FilterWidth();
 
     switch(kindHybridRANSLES){
       case SA_DES: {
@@ -1405,9 +1473,18 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
         1997
         ---*/
 
-        const su2double maxDelta = geometry->nodes->GetMaxLength(iPoint);
+        su2double maxDelta = geometry->nodes->GetMaxLength(iPoint);
+        if (LES_FilterWidth > 0.0){
+          maxDelta = LES_FilterWidth;
+        }
         const su2double distDES = constDES * maxDelta;
         lengthScale = min(distDES,wallDistance);
+        lesSensor = (wallDistance<=distDES) ? 0.0 : 1.0;
+
+        if (config->GetEnforceLES()) {
+          lengthScale = distDES;
+          lesSensor = 1.0;
+        }
 
         break;
       }
@@ -1417,13 +1494,22 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
          Theoretical and Computational Fluid Dynamics - 2006
          ---*/
 
-        const su2double maxDelta = geometry->nodes->GetMaxLength(iPoint);
+        su2double maxDelta = geometry->nodes->GetMaxLength(iPoint);
+        if (LES_FilterWidth > 0.0){
+          maxDelta = LES_FilterWidth;
+        }
 
         const su2double r_d = (kinematicViscosityTurb+kinematicViscosity)/(uijuij*k2*pow(wallDistance, 2.0));
         const su2double f_d = 1.0-tanh(pow(8.0*r_d,3.0));
 
         const su2double distDES = constDES * maxDelta;
         lengthScale = wallDistance-f_d*max(0.0,(wallDistance-distDES));
+        lesSensor = (wallDistance<=distDES) ? 0.0 : f_d;
+
+        if (config->GetEnforceLES()) {
+          lengthScale = distDES;
+          lesSensor = 1.0;
+        }
 
         break;
       }
@@ -1462,8 +1548,17 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
           maxDelta = deltaDDES;
         }
 
+        if (LES_FilterWidth > 0.0){
+          maxDelta = LES_FilterWidth;
+        }
         const su2double distDES = constDES * maxDelta;
         lengthScale = wallDistance-f_d*max(0.0,(wallDistance-distDES));
+        lesSensor = (wallDistance<=distDES) ? 0.0 : f_d;
+
+        if (config->GetEnforceLES()) {
+          lengthScale = distDES;
+          lesSensor = 1.0;
+        }
 
         break;
       }
@@ -1514,17 +1609,286 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
           maxDelta = deltaDDES;
         }
 
+        if (LES_FilterWidth > 0.0){
+          maxDelta = LES_FilterWidth;
+        }
         const su2double distDES = constDES * maxDelta;
         lengthScale = wallDistance-f_d*max(0.0,(wallDistance-distDES));
+        lesSensor = (wallDistance<=distDES) ? 0.0 : f_d;
+
+        if (config->GetEnforceLES()) {
+          lengthScale = distDES;
+          lesSensor = 1.0;
+        }
 
         break;
       }
     }
 
     nodes->SetDES_LengthScale(iPoint, lengthScale);
+    nodes->SetLES_Mode(iPoint, lesSensor);
+
+  }
+  END_SU2_OMP_FOR
+}
+
+void CTurbSASolver::SetLangevinSourceTerms(CConfig *config, CGeometry* geometry) {
+
+  unsigned long timeIter = config->GetTimeIter();
+  const su2double threshold = config->GetStochFdThreshold();
+  const su2double dummySource = 1e3;
+  bool stochBackscatterInBox = config->GetStochBackscatterInBox();
+  bool insideBox = true;
+
+  SU2_OMP_FOR_DYN(omp_chunk_size)
+  for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++){
+    unsigned long iPointGlobal = geometry->nodes->GetGlobalIndex(iPoint);
+    for (unsigned short iDim = 0; iDim < nDim; iDim++){
+      su2double lesSensor = nodes->GetLES_Mode(iPoint);
+      const auto coord = geometry->nodes->GetCoord(iPoint);
+      if (stochBackscatterInBox) {
+        auto sbsBoxBounds = config->GetStochBackscatterBoxBounds();
+        bool insideBoxX = (coord[0]>=sbsBoxBounds[0] && coord[0]<=sbsBoxBounds[1]);
+        bool insideBoxY = (coord[1]>=sbsBoxBounds[2] && coord[1]<=sbsBoxBounds[3]);
+        bool insideBoxZ = (coord[2]>=sbsBoxBounds[4] && coord[2]<=sbsBoxBounds[5]);
+        insideBox = (insideBoxX && insideBoxY && insideBoxZ);
+      }
+      if (lesSensor>threshold && insideBox) {
+        su2double rnd = RandomToolbox::GetNormal(iPointGlobal, iDim, timeIter);
+        nodes->SetLangevinSourceTermsOld(iPoint, iDim, rnd);
+        nodes->SetLangevinSourceTerms(iPoint, iDim, rnd);
+      } else {
+        nodes->SetLangevinSourceTermsOld(iPoint, iDim, dummySource);
+        nodes->SetLangevinSourceTerms(iPoint, iDim, 0.0);
+      }
+    }
+  }
+  END_SU2_OMP_FOR
 
+  SU2_OMP_FOR_DYN(omp_chunk_size)
+  for (unsigned short iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+    for (unsigned long iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) {
+      unsigned long iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
+      if (config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE) {
+        for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+          nodes->SetLangevinSourceTermsOld(iPoint, iDim, dummySource);
+          nodes->SetLangevinSourceTerms(iPoint, iDim, 0.0);
+        }
+      }
+    }
   }
   END_SU2_OMP_FOR
+
+}
+
+void CTurbSASolver::SmoothLangevinSourceTerms(CConfig* config, CGeometry* geometry) {
+
+  const su2double LES_FilterWidth = config->GetLES_FilterWidth();
+  const su2double constDES = config->GetConst_DES();
+  const su2double cDelta = config->GetSBS_Cdelta();
+  const unsigned short maxIter = config->GetSBS_maxIterSmooth();
+  const su2double tol = -5.0;
+  const su2double sourceLim = 3.0;
+  const su2double omega = 0.8;
+  unsigned long timeIter = config->GetTimeIter();
+  unsigned long restartIter = config->GetRestart_Iter();
+
+  /*--- Start SOR algorithm for the Laplacian smoothing. ---*/
+
+  for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+  
+    for (unsigned short iter = 0; iter < maxIter; iter++) {
+      
+      /*--- MPI communication. ---*/
+
+      InitiateComms(geometry, config, MPI_QUANTITIES::STOCH_SOURCE_LANG);
+      CompleteComms(geometry, config, MPI_QUANTITIES::STOCH_SOURCE_LANG);
+
+      su2double localResNorm = 0.0;
+      unsigned long local_nPointLES = 0;
+
+      SU2_OMP_FOR_DYN(omp_chunk_size)
+      for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
+        su2double source_i_old = nodes->GetLangevinSourceTermsOld(iPoint, iDim);
+        if (source_i_old > 3.0*sourceLim) continue;
+        local_nPointLES += 1;
+        const su2double DES_LengthScale = nodes->GetDES_LengthScale(iPoint);
+        su2double maxDelta = DES_LengthScale / constDES;
+        if (LES_FilterWidth > 0.0) maxDelta = LES_FilterWidth;
+        su2double b = sqrt(cDelta) * maxDelta;
+        su2double b2 = b * b;
+        su2double volume_iPoint = geometry->nodes->GetVolume(iPoint);
+        su2double source_i = nodes->GetLangevinSourceTerms(iPoint, iDim);
+        auto coord_i = geometry->nodes->GetCoord(iPoint);
+
+        /*--- Assemble system matrix. ---*/
+
+        su2double diag = 1.0;
+        su2double sum = 0.0;
+        for (unsigned short iNode = 0; iNode < geometry->nodes->GetnPoint(iPoint); iNode++) {
+          auto jPoint = geometry->nodes->GetPoint(iPoint, iNode);
+          auto coord_j = geometry->nodes->GetCoord(jPoint);
+          auto iEdge = geometry->nodes->GetEdge(iPoint, iNode);
+          auto* normal = geometry->edges->GetNormal(iEdge);
+          su2double area = GeometryToolbox::Norm(nDim, normal);
+          su2double dx_ij = GeometryToolbox::Distance(nDim, coord_i, coord_j);
+          su2double source_j = nodes->GetLangevinSourceTerms(jPoint, iDim);
+          su2double a_ij = area/volume_iPoint * b2/dx_ij;
+          diag += a_ij;
+          sum += a_ij * source_j;
+        }
+
+        /*--- Update the solution. ---*/
+
+        su2double source_tmp = (source_i_old + sum) / diag;
+        localResNorm += pow(omega * (source_tmp - source_i), 2);
+        source_i = (1.0-omega)*source_i + omega*source_tmp;
+        nodes->SetLangevinSourceTerms(iPoint, iDim, source_i);
+
+      }
+      END_SU2_OMP_FOR
+
+      /*--- Stop integration if residual drops below tolerance. ---*/
+
+      su2double globalResNorm = 0.0;
+      unsigned long global_nPointLES = 0;
+      SU2_MPI::Allreduce(&local_nPointLES, &global_nPointLES, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
+      SU2_MPI::Allreduce(&localResNorm, &globalResNorm, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+      globalResNorm = (global_nPointLES==0) ? 0.0 : sqrt(globalResNorm / global_nPointLES);
+
+      if (rank == MASTER_NODE) {
+        if (iter == 0) {
+          cout << endl
+               << "Residual of Laplacian smoothing along dimension " << iDim+1 << "." << endl
+               << "---------------------------------" << endl
+               << "   Iter       RMS Residual" << endl
+               << "---------------------------------" << endl;
+        }
+        if (iter%10 == 0) {
+          cout << "  " 
+               << std::setw(5) << iter
+               << "       "
+               << std::setw(12) << std::fixed << std::setprecision(6) << log10(globalResNorm)
+               << endl;
+        }
+      }
+
+      if (log10(globalResNorm) < tol || iter == maxIter-1) {
+
+        if (rank == MASTER_NODE) {
+          cout << "  " 
+               << std::setw(5) << iter
+               << "       "
+               << std::setw(12) << std::fixed << ::setprecision(6) << log10(globalResNorm)
+               << endl;
+          cout << "---------------------------------" << endl;
+        }
+        
+        /*--- Scale source terms for variance preservation. ---*/
+
+        su2double var_check_old = 0.0;
+        su2double mean_check_old = 0.0;
+        su2double var_check_new = 0.0;
+        su2double mean_check_new = 0.0;
+        su2double var_check_notSmoothed = 0.0;
+        su2double mean_check_notSmoothed = 0.0;
+        for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
+          su2double source_notSmoothed = nodes->GetLangevinSourceTermsOld(iPoint, iDim);
+          if (source_notSmoothed > 3.0*sourceLim) continue;
+          su2double source = nodes->GetLangevinSourceTerms(iPoint, iDim);
+          mean_check_old += source;
+          var_check_old += source * source;
+          mean_check_notSmoothed += source_notSmoothed;
+          var_check_notSmoothed += source_notSmoothed * source_notSmoothed;
+          su2double integral = 0.0;
+          if (timeIter==restartIter) {
+            const su2double DES_LengthScale = nodes->GetDES_LengthScale(iPoint);
+            su2double maxDelta = DES_LengthScale / constDES;
+            if (LES_FilterWidth > 0.0) maxDelta = LES_FilterWidth;
+            su2double b = sqrt(cDelta) * maxDelta;
+            su2double b2 = b * b;
+            auto coord_i = geometry->nodes->GetCoord(iPoint);
+            su2double max_dx = 0.0;
+            su2double max_dy = 0.0;
+            su2double max_dz = 0.0;
+            unsigned short nNeigh = geometry->nodes->GetnPoint(iPoint);
+            for (unsigned short iNode = 0; iNode < nNeigh; iNode++) {
+              auto jPoint = geometry->nodes->GetPoint(iPoint, iNode);
+              auto coord_j = geometry->nodes->GetCoord(jPoint);
+              su2double dx = fabs(coord_i[0]-coord_j[0]);
+              su2double dy = fabs(coord_i[1]-coord_j[1]);
+              su2double dz = 0.0;
+              if (nDim == 3) dz = fabs(coord_i[2]-coord_j[2]);
+              if (dx > max_dx) max_dx = dx;
+              if (dy > max_dy) max_dy = dy;
+              if (dz > max_dz) max_dz = dz;
+            }
+            su2double dx2 = max_dx * max_dx;
+            su2double dy2 = max_dy * max_dy;
+            su2double dz2 = max_dz * max_dz;
+            su2double bx = b2 / dx2;
+            su2double by = b2 / dy2;
+            su2double bz = 0.0;
+            if (nDim == 3) bz = b2 / dz2;
+            integral = RandomToolbox::GetBesselIntegral(bx, by, bz);
+            nodes->SetBesselIntegral(iPoint, integral);
+          } else {
+            integral = nodes->GetBesselIntegral(iPoint);
+          }
+          su2double scaleFactor = 1.0 / sqrt(max(integral, 1e-10));
+          source *= scaleFactor;
+          source = min(max(source, -sourceLim), sourceLim);
+          mean_check_new += source;
+          var_check_new += source * source;
+          nodes->SetLangevinSourceTerms(iPoint, iDim, source);
+        }
+        su2double mean_check_new_G = 0.0;
+        SU2_MPI::Allreduce(&mean_check_new, &mean_check_new_G, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+        mean_check_new_G /= global_nPointLES;
+        for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
+          su2double source_notSmoothed = nodes->GetLangevinSourceTermsOld(iPoint, iDim);
+          su2double source = nodes->GetLangevinSourceTerms(iPoint, iDim);
+          if (source_notSmoothed > 3.0*sourceLim) continue;
+          source -= mean_check_new_G;
+          nodes->SetLangevinSourceTerms(iPoint, iDim, source);
+        }
+        if (config->GetStochSourceDiagnostics()) {
+          su2double mean_check_old_G = 0.0;
+          su2double mean_check_notSmoothed_G = 0.0;
+          su2double var_check_old_G = 0.0;
+          su2double var_check_new_G = 0.0;
+          su2double var_check_notSmoothed_G = 0.0;
+          SU2_MPI::Allreduce(&mean_check_old, &mean_check_old_G, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+          SU2_MPI::Allreduce(&mean_check_notSmoothed, &mean_check_notSmoothed_G, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+          SU2_MPI::Allreduce(&var_check_old, &var_check_old_G, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+          SU2_MPI::Allreduce(&var_check_new, &var_check_new_G, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+          SU2_MPI::Allreduce(&var_check_notSmoothed, &var_check_notSmoothed_G, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+          mean_check_old_G /= global_nPointLES;
+          var_check_old_G /= global_nPointLES;
+          var_check_old_G -= mean_check_old_G * mean_check_old_G;
+          var_check_new_G /= global_nPointLES;
+          var_check_new_G -= mean_check_new_G * mean_check_new_G;
+          mean_check_notSmoothed_G /= global_nPointLES;
+          var_check_notSmoothed_G /= global_nPointLES;
+          var_check_notSmoothed_G -= mean_check_notSmoothed_G * mean_check_notSmoothed_G;
+          if (rank == MASTER_NODE) {
+            cout << "Mean of stochastic source term in Langevin equations: " << endl;
+            cout << "   Uncorrelated            --> " << mean_check_notSmoothed_G << endl;
+            cout << "   Smoothed before scaling --> " << mean_check_old_G << endl;
+            cout << "   Smoothed after scaling  --> " << mean_check_new_G << " (subtracted from stochastic field to guarantee zero mean)" << endl;
+            cout << "Variance of stochastic source term in Langevin equations: " << endl;
+            cout << "   Uncorrelated            --> " << var_check_notSmoothed_G << endl;
+            cout << "   Smoothed before scaling --> " << var_check_old_G << endl;
+            cout << "   Smoothed after scaling  --> " << var_check_new_G << endl;
+            cout << endl;
+          }
+        }
+        break;
+
+      }
+    }
+  }
+
 }
 
 void CTurbSASolver::SetInletAtVertex(const su2double *val_inlet,
@@ -1561,7 +1925,8 @@ void CTurbSASolver::ComputeUnderRelaxationFactor(const CConfig *config) {
    SA_NEG model is more robust due to allowing for negative nu_tilde,
    so the under-relaxation is not applied to that variant. */
 
-  if (config->GetSAParsedOptions().version == SA_OPTIONS::NEG) return;
+  if (config->GetSAParsedOptions().version == SA_OPTIONS::NEG ||
+      config->GetStochastic_Backscatter()) return;
 
   /* Loop over the solution update given by relaxing the linear
    system for this nonlinear iteration. */
diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp
index f5201c92cb58..fd2cb97c60bc 100644
--- a/SU2_CFD/src/variables/CIncNSVariable.cpp
+++ b/SU2_CFD/src/variables/CIncNSVariable.cpp
@@ -37,6 +37,7 @@ CIncNSVariable::CIncNSVariable(su2double pressure, const su2double *velocity, su
   StrainMag.resize(nPoint);
   Tau_Wall.resize(nPoint) = su2double(-1.0);
   DES_LengthScale.resize(nPoint) = su2double(0.0);
+  lesMode.resize(nPoint) = su2double(0.0);
   Max_Lambda_Visc.resize(nPoint);
   /*--- Allocate memory for the AuxVar and its gradient. See e.g. CIncEulerSolver::Source_Residual:
    * Axisymmetric: total-viscosity * y-vel / y-coord
diff --git a/SU2_CFD/src/variables/CNEMONSVariable.cpp b/SU2_CFD/src/variables/CNEMONSVariable.cpp
index 37b020c0ba69..3ae287e9d4fc 100644
--- a/SU2_CFD/src/variables/CNEMONSVariable.cpp
+++ b/SU2_CFD/src/variables/CNEMONSVariable.cpp
@@ -70,6 +70,7 @@ CNEMONSVariable::CNEMONSVariable(su2double val_pressure,
   StrainMag.resize(nPoint)       = su2double(0.0);
   Tau_Wall.resize(nPoint)        = su2double(-1.0);
   DES_LengthScale.resize(nPoint) = su2double(0.0);
+  lesMode.resize(nPoint)        = su2double(0.0);
   Roe_Dissipation.resize(nPoint) = su2double(0.0);
   Vortex_Tilting.resize(nPoint)  = su2double(0.0);
   Max_Lambda_Visc.resize(nPoint) = su2double(0.0);
diff --git a/SU2_CFD/src/variables/CNSVariable.cpp b/SU2_CFD/src/variables/CNSVariable.cpp
index 365c5b5d870c..951c1c0cff8b 100644
--- a/SU2_CFD/src/variables/CNSVariable.cpp
+++ b/SU2_CFD/src/variables/CNSVariable.cpp
@@ -39,6 +39,7 @@ CNSVariable::CNSVariable(su2double density, const su2double *velocity, su2double
   StrainMag.resize(nPoint) = su2double(0.0);
   Tau_Wall.resize(nPoint) = su2double(-1.0);
   DES_LengthScale.resize(nPoint) = su2double(0.0);
+  lesMode.resize(nPoint) = su2double(0.0);
   Roe_Dissipation.resize(nPoint) = su2double(0.0);
   Vortex_Tilting.resize(nPoint) = su2double(0.0);
   Max_Lambda_Visc.resize(nPoint) = su2double(0.0);
diff --git a/SU2_CFD/src/variables/CTurbSAVariable.cpp b/SU2_CFD/src/variables/CTurbSAVariable.cpp
index de9c6d6bb481..8f107936725e 100644
--- a/SU2_CFD/src/variables/CTurbSAVariable.cpp
+++ b/SU2_CFD/src/variables/CTurbSAVariable.cpp
@@ -28,12 +28,22 @@
 
 #include "../../include/variables/CTurbSAVariable.hpp"
 
-
 CTurbSAVariable::CTurbSAVariable(su2double val_nu_tilde, su2double val_muT, unsigned long npoint,
                                  unsigned long ndim, unsigned long nvar, CConfig *config) :
                  CTurbVariable(npoint, ndim, nvar, config) {
 
-  Solution_Old = Solution = val_nu_tilde;
+  /*--- Initialize solution (check if the Stochastic Backscatter Model is active) ---*/
+  bool backscatter = config->GetStochastic_Backscatter();
+  if (!backscatter) {
+    Solution_Old = Solution = val_nu_tilde;
+  } else {
+    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+      Solution_Old(iPoint, 0) = Solution(iPoint, 0) = val_nu_tilde;
+      for (unsigned long iVar = 1; iVar < nVar; iVar++) {
+        Solution_Old(iPoint, iVar) = Solution(iPoint, iVar) = 0.0;
+      }
+    }
+  }
 
   muT.resize(nPoint) = val_muT;
 
@@ -47,7 +57,12 @@ CTurbSAVariable::CTurbSAVariable(su2double val_nu_tilde, su2double val_muT, unsi
   }
 
   DES_LengthScale.resize(nPoint) = su2double(0.0);
+  lesMode.resize(nPoint) = su2double(0.0);
   Vortex_Tilting.resize(nPoint);
+  stochSource.resize(nPoint, nDim) = su2double(0.0);
+  stochSourceOld.resize(nPoint, nDim) = su2double(0.0);
+  besselIntegral.resize(nPoint);
+
 }
 
 void CTurbSAVariable::SetVortex_Tilting(unsigned long iPoint, CMatrixView<const su2double> PrimGrad_Flow,
diff --git a/TestCases/backscatter/DIHT/backscatter_DIHT.cfg b/TestCases/backscatter/DIHT/backscatter_DIHT.cfg
new file mode 100644
index 000000000000..6fd2aa724855
--- /dev/null
+++ b/TestCases/backscatter/DIHT/backscatter_DIHT.cfg
@@ -0,0 +1,107 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: Decaying Isotropic Homogeneous Turbulence (DIHT)           % 
+% Author: Angelo Passariello                                                   %
+% Institution: University of Naples Federico II                                %
+% Date: 2025.10.01                                                             %
+% File Version 8.3.0 "Harrier"                                                 %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+
+SOLVER= RANS
+KIND_TURB_MODEL= SA
+HYBRID_RANSLES= SA_DES 
+MATH_PROBLEM= DIRECT
+RESTART_SOL= NO
+
+% ------------- STOCHASTIC BACKSCATTER MODEL PARAMETERS -----------------------%
+
+STOCHASTIC_BACKSCATTER= YES
+SBS_CTAU= 0.001
+
+% ------------- DECAYING ISOTROPIC HOMOGENEOUS TURBULENCE ---------------------%
+
+DIHT= YES
+DIHT_DOMAIN_LENGTH= (0.5588, 0.5588, 0.5588)
+DIHT_NPOINT= (64, 64, 64) 
+DIHT_NUM_MODES= 800
+
+% ----------- COMPRESSIBLE AND INCOMPRESSIBLE FREE-STREAM DEFINITION ----------%
+
+MACH_NUMBER= 0.001
+FREESTREAM_TEMPERATURE= 300.0
+REYNOLDS_NUMBER= 10129
+REYNOLDS_LENGTH= 0.5588
+FREESTREAM_NU_FACTOR= 1.0
+
+% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
+
+REF_LENGTH= 1.0
+REF_AREA= 1.0
+
+% ------------------------- UNSTEADY SIMULATION -------------------------------%
+
+TIME_DOMAIN= NO
+%TIME_MARCHING= DUAL_TIME_STEPPING-2ND_ORDER
+%TIME_STEP= 0.001
+%MAX_TIME= 20
+%UNST_CFL_NUMBER= 0.0
+INNER_ITER= 0  
+
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+
+MARKER_PERIODIC= (xmin, xmax, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5588, 0.0, 0.0, ymin, ymax, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5588, 0.0, zmin, zmax, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5588)
+
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+
+NUM_METHOD_GRAD= GREEN_GAUSS
+CFL_NUMBER= 10.0
+CFL_ADAPT= NO
+CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 )
+RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 )
+TIME_ITER= 1
+LINEAR_SOLVER_PREC= ILU
+LOW_MACH_PREC= YES
+%LOW_MACH_CORR= YES
+
+% ----------------------- SLOPE LIMITER DEFINITION ----------------------------%
+
+MUSCL_FLOW= YES
+SLOPE_LIMITER_FLOW= NONE
+MUSCL_TURB= NO
+VENKAT_LIMITER_COEFF= 0.05
+
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+CONV_NUM_METHOD_FLOW= SLAU2
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+
+% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------%
+%
+CONV_NUM_METHOD_TURB= SCALAR_UPWIND
+TIME_DISCRE_TURB= EULER_IMPLICIT
+
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+CONV_RESIDUAL_MINVAL= -15
+CONV_STARTITER= 0
+CONV_CAUCHY_ELEMS= 100
+CONV_CAUCHY_EPS= 1E-6
+
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+MESH_FILENAME= cube_64.su2 
+MESH_FORMAT= SU2
+TABULAR_FORMAT= CSV
+CONV_FILENAME= history
+RESTART_FILENAME= restart_flow
+VOLUME_FILENAME= flow
+VOLUME_OUTPUT= DDES, BACKSCATTER
+SCREEN_OUTPUT= (TIME_ITER, INNER_ITER, RMS_DENSITY, RMS_MOMENTUM-X, RMS_MOMENTUM-Y, RMS_MOMENTUM-Z, RMS_NU_TILDE, RMS_STOCH_VAR_X, RMS_STOCH_VAR_Y, RMS_STOCH_VAR_Z)
+HISTORY_OUTPUT= (TIME_ITER, RMS_DENSITY, RMS_MOMENTUM-X, RMS_MOMENTUM-Y, RMS_MOMENTUM-Z, RMS_NU_TILDE)
+HISTORY_WRT_FREQ_INNER= 1000000
+OUTPUT_WRT_FREQ= 1000000
+OUTPUT_FILES= (RESTART, PARAVIEW)
diff --git a/TestCases/ddes/flatplate/ddes_flatplate.cfg b/TestCases/ddes/flatplate/ddes_flatplate.cfg
index a65727ab04d1..723c11fcb574 100644
--- a/TestCases/ddes/flatplate/ddes_flatplate.cfg
+++ b/TestCases/ddes/flatplate/ddes_flatplate.cfg
@@ -13,9 +13,11 @@
 %
 SOLVER= RANS
 KIND_TURB_MODEL= SA
-HYBRID_RANSLES= SA_EDDES
+HYBRID_RANSLES= SA_DES
+STOCHASTIC_BACKSCATTER= YES
 MATH_PROBLEM= DIRECT
 RESTART_SOL= NO
+RESTART_ITER= 100
 
 % ----------- COMPRESSIBLE AND INCOMPRESSIBLE FREE-STREAM DEFINITION ----------%
 %
@@ -41,7 +43,7 @@ TIME_MARCHING= DUAL_TIME_STEPPING-2ND_ORDER
 %
 % U_inf = 69.4448 - dt*=0.02 - dt=0.000288
 TIME_STEP= 0.000288
-MAX_TIME= 20.0
+MAX_TIME= 5
 UNST_CFL_NUMBER= 0.0
 INNER_ITER= 20
 
@@ -61,7 +63,7 @@ CFL_NUMBER= 10.0
 CFL_ADAPT= NO
 CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 )
 RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 )
-TIME_ITER= 10
+TIME_ITER= 500
 
 % ----------------------- SLOPE LIMITER DEFINITION ----------------------------%
 %
@@ -103,5 +105,5 @@ VOLUME_ADJ_FILENAME= adjoint
 GRAD_OBJFUNC_FILENAME= of_grad
 SURFACE_FILENAME= surface_flow
 SURFACE_ADJ_FILENAME= surface_adjoint
-OUTPUT_WRT_FREQ= 1000
+OUTPUT_WRT_FREQ= 1000000
 SCREEN_OUTPUT= (TIME_ITER, INNER_ITER, RMS_DENSITY, RMS_NU_TILDE, LIFT, DRAG, TOTAL_HEATFLUX)
diff --git a/config_template.cfg b/config_template.cfg
index 99eb260e61ed..468ccadcc127 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -180,7 +180,43 @@ HYBRID_RANSLES= SA_DDES
 %
 % DES Constant (0.65)
 DES_CONST= 0.65
-
+%
+% User-specified LES filter width (if negative, it is computed based on the local cell size, default -1.0)
+LES_FILTER_WIDTH= -1.0
+%
+% Stochastic Backscatter Model (NO, YES)
+STOCHASTIC_BACKSCATTER= NO
+%
+% Backscatter lengthscale coefficient (0.1)
+SBS_LENGTHSCALE_COEFF= 0.1
+%
+% Maximum number of smoothing iterations for SBS model (100)
+SBS_MAX_ITER_SMOOTH= 100
+%
+% Backscatter timescale coefficient (0.05)
+SBS_TIMESCALE_COEFF= 0.05
+%
+% Backscatter intensity coefficient (1.0)
+SBS_INTENSITY_COEFF= 1.0
+%
+% Relaxation factor for the stochastic source term in the main balance equations (0.0)
+SBS_RELAXATION_FACTOR= 0.0
+%
+% Include stochastic source in turbulence model equation (NO, YES)
+SBS_SOURCE_NU_EQUATION= YES
+%
+% Include diagnostics of the stochastic source term in Langevin equations (NO, YES)
+SBS_SOURCE_DIAGNOSTICS= NO
+%
+% Apply Stochastic Backscatter Model only in a bounded box (NO, YES)
+SBS_IN_BOX= NO
+%
+% Bounds of the box where the Stochastic Backscatter Model is activated (x_min, x_max, y_min, y_max, z_min, z_max)
+SBS_BOX_BOUNDS= (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+%
+% Enforce LES mode in Hybrid RANS-LES Simulations (NO, YES)
+ENFORCE_LES= NO
+%
 % -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------%
 %
 % Mach number (non-dimensional, based on the free-stream values)
@@ -1640,7 +1676,7 @@ SMOOTH_GEOMETRY= 0
 %
 % Convective numerical method (JST, JST_KE, JST_MAT, LAX-FRIEDRICH, ROE, AUSM,
 %                              AUSMPLUSUP, AUSMPLUSUP2, AUSMPLUSM, HLLC, TURKEL_PREC,
-%                              SW, MSW, FDS, SLAU, SLAU2, L2ROE, LMROE)
+%                              SW, MSW, FDS, SLAU, SLAU2, L2ROE, LMROE, LD2)
 CONV_NUM_METHOD_FLOW= ROE
 %
 % Roe Low Dissipation function for Hybrid RANS/LES simulations (FD, NTS, NTS_DUCROS)