diff --git a/src/incflo.H b/src/incflo.H
index 82d38b28..b1ffadef 100644
--- a/src/incflo.H
+++ b/src/incflo.H
@@ -38,7 +38,7 @@ public:
 
     enum struct FluidModel {
         Newtonian, powerlaw, Bingham, HerschelBulkley, deSouzaMendesDutra,
-        DataDrivenMPMD
+        DataDrivenMPMD, Rauter
     };
 
     incflo ();
@@ -551,6 +551,7 @@ private:
     bool m_two_fluid = false;
     bool m_two_fluid_eta_harmonic = true;
     amrex::Real m_ro_0_second = 1.0;
+    amrex::Real m_ro_grain_second = 1.0; // Will be used to calculate Inertial number
     FluidModel m_fluid_model_second;
     amrex::Real m_mu_second = 1.0;
     amrex::Real m_n_0_second = 0.0;
@@ -561,6 +562,10 @@ private:
     amrex::Real m_eta_max_second = amrex::Real(1.0e+20);
     amrex::Real m_diam_second = amrex::Real(0.);
     amrex::Real m_min_conc_second = amrex::Real(0.); // Used to limit expensive particle-sims
+    // These are specifically for Rauter_2021
+    amrex::Real m_mu_1_second = amrex::Real(0.);
+    amrex::Real m_mu_2_second = amrex::Real(0.);
+    amrex::Real m_I_0_second = amrex::Real(0.);
 
     int m_plot_int = -1;
 
diff --git a/src/rheology/incflo_read_rheology_parameters.cpp b/src/rheology/incflo_read_rheology_parameters.cpp
index 816fc7b3..d42ed7d2 100644
--- a/src/rheology/incflo_read_rheology_parameters.cpp
+++ b/src/rheology/incflo_read_rheology_parameters.cpp
@@ -107,6 +107,9 @@ void incflo::ReadRheologyParameters()
         amrex::ParmParse pp_scnd("incflo.second_fluid");
         pp_scnd.query("ro_0", m_ro_0_second);
         AMREX_ALWAYS_ASSERT(m_ro_0_second >= 0.0);
+        // Initially setting ro_grain the same as ro_0
+        m_ro_grain_second = m_ro_0_second;
+        pp_scnd.query("ro_grain",m_ro_grain_second); 
         pp_scnd.query("mu", m_mu_second);
         std::string fluid_model_s_snd = "newtonian";
         pp_scnd.query("fluid_model", fluid_model_s_snd);
@@ -201,6 +204,14 @@ void incflo::ReadRheologyParameters()
             amrex::Print() << "Data-driven model through AMReX-MPMD."<<std::endl;
         }
 #endif
+        else if(fluid_model_s_snd == "rauter")
+        {
+            m_fluid_model_second = FluidModel::Rauter;
+            pp_scnd.get("mu_1", m_mu_1_second);
+            pp_scnd.get("mu_2", m_mu_2_second);
+            pp_scnd.get("I_0", m_I_0_second);
+            amrex::Print() << "Using mu(I) defined Rauter 2021 (Eq. 2.30)"<<std::endl;
+        }
         else
         {
             amrex::Abort("Unknown fluid_model! Choose either newtonian, powerlaw, bingham, hb, smd");
diff --git a/src/rheology/incflo_rheology.cpp b/src/rheology/incflo_rheology.cpp
index 369bde73..be59c4f4 100644
--- a/src/rheology/incflo_rheology.cpp
+++ b/src/rheology/incflo_rheology.cpp
@@ -404,7 +404,7 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/,
            // NOTE: HYDROSTATIC PRESSURE IS INCORRECT IF THERE ARE
            // MULTIPLE BOXES/GRIDS ALONG GRAVITY DIRECTION
 
-           // Inertial Number
+           // Inertial Number = diameter*strainrate*sqrt(rho_grain/p)
            // NOTE: Strain-rate calculated is TWO TIMES the actual value
            // The second component will carry concentration
            MultiFab inertial_num(vel_eta->boxArray(),vel_eta->DistributionMap(),2,nghost);
@@ -419,11 +419,11 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/,
                Array4<Real> const& inrt_num_arr = inertial_num.array(mfi);
                const Real eps = Real(1.0e-20);
                const Real diam_scnd = m_diam_second;
-               const Real ro_0_scnd = m_ro_0_second;
+               const Real ro_scnd = m_ro_grain_second;
                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                {
                     inrt_num_arr(i,j,k,0) =
-                       std::sqrt(ro_0_scnd/(p_static_arr(i,j,k)+eps))*
+                       std::sqrt(ro_scnd/(p_static_arr(i,j,k)+eps))*
                        diam_scnd*Real(0.5)*sr_arr(i,j,k);
                });
            }
@@ -456,6 +456,37 @@ void incflo::compute_nodal_viscosity_at_level (int /*lev*/,
 
            } else
 #endif
+           if (m_fluid_model_second == FluidModel::Rauter) {
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+               for (MFIter mfi(sr_mf,TilingIfNotGPU()); mfi.isValid(); ++mfi)
+               {
+                   Box const& bx = mfi.growntilebox(nghost);
+                   Array4<Real const> const& sr_arr = sr_mf.const_array(mfi);
+                   Array4<Real const> const& p_static_arr = p_static.const_array(mfi);
+                   Array4<Real const> const& inrt_num_arr = inertial_num.const_array(mfi);
+                   Array4<Real> const& vel_eta_snd_arr = vel_eta_second.array(mfi);
+                   const Real mu_1_scnd = m_mu_1_second;
+                   const Real mu_2_scnd = m_mu_2_second;
+                   const Real I_0_scnd = m_I_0_second;
+                   const Real eps = Real(1.0e-20);
+                   // Note: sr_mf contains TWO TIMES strain rate
+                   // Note: Inertial number in Rauter 2021 (Eq. 2.29)
+                   // is 2 times the calculated value (inertial_num)
+                   amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                   {
+                        vel_eta_snd_arr(i,j,k) = (inrt_num_arr(i,j,k,0)*Real(2.0));
+                        vel_eta_snd_arr(i,j,k) /= (I_0_scnd + Real(2.0)*inrt_num_arr(i,j,k,0));
+                        vel_eta_snd_arr(i,j,k) *= (mu_2_scnd-mu_1_scnd);
+                        vel_eta_snd_arr(i,j,k) += mu_1_scnd;
+                        // The above value is stress ratio
+                        vel_eta_snd_arr(i,j,k) *= p_static_arr(i,j,k);
+                        vel_eta_snd_arr(i,j,k) /= (sr_arr(i,j,k)+eps);
+                   });
+               }
+
+           } else
            {
                NonNewtonianViscosity non_newtonian_viscosity;
                non_newtonian_viscosity.fluid_model = m_fluid_model_second;