diff --git a/include/node_kernels/SDRSSTAMSNodeKernel.h b/include/node_kernels/SDRSSTAMSNodeKernel.h
index 7c8ee4258..9b9f0ebb7 100644
--- a/include/node_kernels/SDRSSTAMSNodeKernel.h
+++ b/include/node_kernels/SDRSSTAMSNodeKernel.h
@@ -79,6 +79,7 @@ class SDRSSTAMSNodeKernel : public NGPNodeKernel<SDRSSTAMSNodeKernel>
   NodeKernelTraits::DblType referenceVelocity_;
   NodeKernelTraits::DblType gammaOne_;
   NodeKernelTraits::DblType gammaTwo_;
+  NodeKernelTraits::DblType sdrAmb_;
 };
 
 } // namespace nalu
diff --git a/include/node_kernels/TKESSTAMSNodeKernel.h b/include/node_kernels/TKESSTAMSNodeKernel.h
index 722dca908..2ce1c1e18 100644
--- a/include/node_kernels/TKESSTAMSNodeKernel.h
+++ b/include/node_kernels/TKESSTAMSNodeKernel.h
@@ -62,6 +62,8 @@ class TKESSTAMSNodeKernel : public NGPNodeKernel<TKESSTAMSNodeKernel>
 
   NodeKernelTraits::DblType betaStar_;
   NodeKernelTraits::DblType tkeProdLimitRatio_;
+  NodeKernelTraits::DblType tkeAmb_;
+  NodeKernelTraits::DblType sdrAmb_;
   const int nDim_;
 };
 
diff --git a/src/node_kernels/SDRSSTAMSNodeKernel.C b/src/node_kernels/SDRSSTAMSNodeKernel.C
index c2b25d3db..a54b809f8 100644
--- a/src/node_kernels/SDRSSTAMSNodeKernel.C
+++ b/src/node_kernels/SDRSSTAMSNodeKernel.C
@@ -62,6 +62,7 @@ SDRSSTAMSNodeKernel::setup(Realm& realm)
   gammaOne_ = realm.get_turb_model_constant(TM_gammaOne);
   gammaTwo_ = realm.get_turb_model_constant(TM_gammaTwo);
   tkeProdLimitRatio_ = realm.get_turb_model_constant(TM_tkeProdLimitRatio);
+  sdrAmb_ = realm.get_turb_model_constant(TM_sdrAmb);
 
   lengthScaleLimiter_ = realm.solutionOptions_->lengthScaleLimiter_;
   if (lengthScaleLimiter_) {
@@ -146,7 +147,10 @@ SDRSSTAMSNodeKernel::execute(
 
   const NodeKernelTraits::DblType dualVolume = dualNodalVolume_.get(node, 0);
 
-  rhs(0) += (Pw - Dw + Sw) * dualVolume;
+  // SUST source term
+  const NodeKernelTraits::DblType Dwamb = beta * rho * sdrAmb_ * sdrAmb_;
+
+  rhs(0) += (Pw - Dw + Dwamb + Sw) * dualVolume;
 
   lhs(0, 0) +=
     (2.0 * beta * rho * sdr + stk::math::max(Sw / sdr, 0.0)) * dualVolume;
diff --git a/src/node_kernels/TKESSTAMSNodeKernel.C b/src/node_kernels/TKESSTAMSNodeKernel.C
index 1aac43bc1..910cee463 100644
--- a/src/node_kernels/TKESSTAMSNodeKernel.C
+++ b/src/node_kernels/TKESSTAMSNodeKernel.C
@@ -51,6 +51,8 @@ TKESSTAMSNodeKernel::setup(Realm& realm)
   // Update turbulence model constants
   betaStar_ = realm.get_turb_model_constant(TM_betaStar);
   tkeProdLimitRatio_ = realm.get_turb_model_constant(TM_tkeProdLimitRatio);
+  tkeAmb_ = realm.get_turb_model_constant(TM_tkeAmb);
+  sdrAmb_ = realm.get_turb_model_constant(TM_sdrAmb);
 }
 
 KOKKOS_FUNCTION
@@ -71,7 +73,11 @@ TKESSTAMSNodeKernel::execute(
 
   const NodeKernelTraits::DblType dualVolume = dualNodalVolume_.get(node, 0);
 
-  rhs(0) += (Pk - Dk) * dualVolume;
+  // SUST source term
+  const NodeKernelTraits::DblType Dkamb =
+    betaStar_ * rho_.get(node, 0) * sdrAmb_ * tkeAmb_;
+
+  rhs(0) += (Pk - Dk + Dkamb) * dualVolume;
 
   lhs(0, 0) += tkeFac * dualVolume;
 }