diff --git a/doc/src/fix_rheo.rst b/doc/src/fix_rheo.rst
index 5eae617419..25e171a1b9 100644
--- a/doc/src/fix_rheo.rst
+++ b/doc/src/fix_rheo.rst
@@ -20,9 +20,10 @@ Syntax
 
        *shift* values = none, turns on velocity shifting
        *thermal* values = none, turns on thermal evolution
-       *surface/detection* values = *sdstyle* *limit*
+       *surface/detection* values = *sdstyle* *limit* *limit/splash*
          *sdstyle* = *coordination* or *divergence*
          *limit* = threshold for surface particles (unitless)
+         *limit/splash* = threshold for splash particles (unitless)
        *interface/reconstruct* values = none, reconstructs interfaces with solid particles
        *rho/sum* values = none, uses the kernel to compute the density of particles
        *density* values = *rho0* (density)
diff --git a/src/RHEO/compute_rheo_surface.cpp b/src/RHEO/compute_rheo_surface.cpp
index 4b0745abd3..230bca8219 100644
--- a/src/RHEO/compute_rheo_surface.cpp
+++ b/src/RHEO/compute_rheo_surface.cpp
@@ -82,6 +82,8 @@ void ComputeRHEOSurface::init()
   threshold_style = fix_rheo->surface_style;
   threshold_divr = fix_rheo->divr_surface;
   threshold_z = fix_rheo->zmin_surface;
+  threshold_splash = fix_rheo->zmin_splash;
+  interface_flag = fix_rheo->interface_flag;
 
   cutsq = cut * cut;
 
@@ -202,13 +204,15 @@ void ComputeRHEOSurface::compute_peratom()
         rhoj = rho[j];
 
         // Add corrections for walls
-        if (fluidi && (!fluidj)) {
-          rhoj = compute_interface->correct_rho(j, i);
-        } else if ((!fluidi) && fluidj) {
-          rhoi = compute_interface->correct_rho(i, j);
-        } else if ((!fluidi) && (!fluidj)) {
-          rhoi = rho0;
-          rhoj = rho0;
+        if (interface_flag) {
+          if (fluidi && (!fluidj)) {
+            rhoj = compute_interface->correct_rho(j, i);
+          } else if ((!fluidi) && fluidj) {
+            rhoi = compute_interface->correct_rho(i, j);
+          } else if ((!fluidi) && (!fluidj)) {
+            rhoi = rho0;
+            rhoj = rho0;
+          }
         }
 
         Voli = mass[itype] / rhoi;
@@ -262,7 +266,7 @@ void ComputeRHEOSurface::compute_peratom()
         if (divr[i] < threshold_divr) {
           status[i] |= STATUS_SURFACE;
           rsurface[i] = 0.0;
-          if (coordination[i] < threshold_z)
+          if (coordination[i] < threshold_splash)
             status[i] |= STATUS_SPLASH;
         }
       }
@@ -272,10 +276,10 @@ void ComputeRHEOSurface::compute_peratom()
       if (mask[i] & groupbit) {
         status[i] |= STATUS_BULK;
         rsurface[i] = cut;
-        if (coordination[i] < threshold_divr) {
+        if (coordination[i] < threshold_z) {
           status[i] |= STATUS_SURFACE;
           rsurface[i] = 0.0;
-          if (coordination[i] < threshold_z)
+          if (coordination[i] < threshold_splash)
             status[i] |= STATUS_SPLASH;
         }
       }
diff --git a/src/RHEO/compute_rheo_surface.h b/src/RHEO/compute_rheo_surface.h
index 220f8beb6d..1d1cdba3fa 100644
--- a/src/RHEO/compute_rheo_surface.h
+++ b/src/RHEO/compute_rheo_surface.h
@@ -41,7 +41,7 @@ class ComputeRHEOSurface : public Compute {
 
  private:
   double cut, cutsq, rho0, threshold_divr;
-  int surface_style, nmax_store, threshold_z;
+  int surface_style, nmax_store, threshold_z, threshold_splash, interface_flag;
   double **B, **gradC;
   int threshold_style, comm_stage;
 
diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp
index 0420a0f23f..339d7f5ac8 100644
--- a/src/RHEO/fix_rheo.cpp
+++ b/src/RHEO/fix_rheo.cpp
@@ -97,14 +97,16 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
       if (strcmp(arg[iarg + 1], "coordination")) {
         surface_style = COORDINATION;
         zmin_surface = utils::inumeric(FLERR,arg[iarg + 2],false,lmp);
+        zmin_splash = utils::inumeric(FLERR,arg[iarg + 3],false,lmp);
       } else if (strcmp(arg[iarg + 1], "divergence")) {
         surface_style = DIVR;
         divr_surface = utils::numeric(FLERR,arg[iarg + 2],false,lmp);
+        zmin_splash = utils::inumeric(FLERR,arg[iarg + 3],false,lmp);
       } else {
         error->all(FLERR,"Illegal surface/detection option in fix rheo, {}", arg[iarg + 1]);
       }
 
-      iarg += 2;
+      iarg += 3;
     } else if (strcmp(arg[iarg],"interface/reconstruct") == 0) {
       interface_flag = 1;
     } else if (strcmp(arg[iarg],"rho/sum") == 0) {
diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h
index 0dbc8db78b..96f4b6f502 100644
--- a/src/RHEO/fix_rheo.h
+++ b/src/RHEO/fix_rheo.h
@@ -40,7 +40,7 @@ class FixRHEO : public Fix {
 
   // Model parameters
   double h, cut, rho0, csq;
-  int zmin_kernel, zmin_surface;
+  int zmin_kernel, zmin_surface, zmin_splash;
   int kernel_style, surface_style;
   double divr_surface;