diff --git a/src/incflo.H b/src/incflo.H
index a3de11ce..c50bdbea 100644
--- a/src/incflo.H
+++ b/src/incflo.H
@@ -593,6 +593,8 @@ private:
     amrex::Real m_mu_p_eps_second = amrex::Real(1.0e-4);
     // epsilon value to strain rate magnitude while calculating dynamic viscosity
     amrex::Real m_mu_sr_eps_second = amrex::Real(1.0e-3);
+    // Variable for surface pressure in mu(I) model
+    amrex::Real m_mu_p_surf_second = amrex::Real(0.);
     // These are related to level set method
     // Coefficient d that is used in level set
     amrex::Real m_level_set_d = amrex::Real(0.1);
diff --git a/src/rheology/incflo_rheology.cpp b/src/rheology/incflo_rheology.cpp
index 4a97b535..477809b5 100644
--- a/src/rheology/incflo_rheology.cpp
+++ b/src/rheology/incflo_rheology.cpp
@@ -227,7 +227,8 @@ void incflo::compute_nodal_viscosity_at_level (int lev,
            compute_nodal_strainrate_at_level(lev,&sr_mf,vel,lev_geom,time,nghost);
            // nodal MultiFab for hydrostatic pressure
            MultiFab p_static(vel_eta->boxArray(),vel_eta->DistributionMap(),1,nghost);
-           compute_nodal_hydrostatic_pressure_at_level(lev,&p_static,rho,Real(0.0),
+           compute_nodal_hydrostatic_pressure_at_level(lev,&p_static,rho,
+                                                       m_mu_p_surf_second,
                                                        lev_geom,nghost);
 
 #ifdef USE_AMREX_MPMD
diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp
index adf94af7..c1944b74 100644
--- a/src/utilities/io.cpp
+++ b/src/utilities/io.cpp
@@ -729,7 +729,7 @@ void incflo::WritePlotFile()
 
             compute_nodal_hydrostatic_pressure_at_level(lev,&p_static,
                                 &m_leveldata[lev]->density,
-                                Real(0.0),
+                                m_mu_p_surf_second,
                                 Geom(lev),0);
 
             MultiFab strainrate(amrex::convert(mf[lev].boxArray(),
@@ -779,7 +779,7 @@ void incflo::WritePlotFile()
 
             compute_nodal_hydrostatic_pressure_at_level(lev,&p_static,
                                 &m_leveldata[lev]->density,
-                                Real(0.0),
+                                m_mu_p_surf_second,
                                 Geom(lev),0);
 
             MultiFab strainrate(amrex::convert(mf[lev].boxArray(),
@@ -818,7 +818,7 @@ void incflo::WritePlotFile()
 
             compute_nodal_hydrostatic_pressure_at_level(lev,&p_static,
                                 &m_leveldata[lev]->density,
-                                Real(0.0),
+                                m_mu_p_surf_second,
                                 Geom(lev),0);
             amrex::average_node_to_cellcenter(mf[lev],icomp,p_static,0,1);
         }