diff --git a/include/nash_cascade.h b/include/nash_cascade.h
index 2a5ad4b..7b809f1 100644
--- a/include/nash_cascade.h
+++ b/include/nash_cascade.h
@@ -11,18 +11,17 @@ struct nash_cascade_parameters {
   double K_nash;             // Fraction of storage per hour that moves from one reservoir to the next (time constant); [1/hour]
   int    nsubsteps;          // the number of substeps that each dt is divided into
   double *nash_storage;      // storage array nash cascade reservoirs; [m]
-  double retention_depth;    // parameter that represents retention depth process that is equivalent [m]
-  // to that used in WRF-Hydro
+  double retention_depth;    /* parameter that represents retention depth process that is equivalent [m]
+				to that used in WRF-Hydro */
   double runon_infiltration; // infiltration losses from surface runoff water to soil (or riparian groundwater) [m/hr]
   int    is_riparian_gw;     // flag to turn on/off riparian groundwater (currently used in LASAM only)
   double K_infiltration;     // Fraction of storage per hour that moves from reservoirs to soil (time constant); [1/hour]
-  double Kinf_c0;            // constant c0 in K_infiltration = c0 + c1 * S, where S is the reservoir storage
-  double Kinf_c1;            // constant c1 in K_infiltration = c0 + c1 * S, where S is the reservoir storage
 };
 
 extern double nash_cascade(double flux_lat_m,int num_lateral_flow_nash_reservoirs,
                            double K_nash,double *nash_storage_arr);
 
-double nash_cascade_surface_runoff(double runoff_m, struct nash_cascade_parameters *nash_params);
+double nash_cascade_surface_runoff(double runoff_m, double soil_storage_deficit_m,
+				   struct nash_cascade_parameters *nash_params);
 
 #endif
diff --git a/src/bmi_cfe.c b/src/bmi_cfe.c
index 9545a95..191161c 100644
--- a/src/bmi_cfe.c
+++ b/src/bmi_cfe.c
@@ -18,14 +18,14 @@
 #define STATE_VAR_NAME_COUNT 94   // must match var_info array size
 
 
-#define PARAM_VAR_NAME_COUNT 19
+#define PARAM_VAR_NAME_COUNT 18
 // NOTE: If you update the params, also update the unit test in ../test/main_unit_test_bmi.c
 static const char *param_var_names[PARAM_VAR_NAME_COUNT] = {
     "maxsmc", "satdk", "slope", "b", "Klf",
     "Kn", "Cgw", "expon", "max_gw_storage",
     "satpsi","wltsmc","alpha_fc","refkdt",
     "a_Xinanjiang_inflection_point_parameter","b_Xinanjiang_shape_parameter","x_Xinanjiang_shape_parameter",
-    "Kinf_c0", "Kinf_c1",
+    "K_infiltration",
     "N_nash"
 };
 
@@ -34,7 +34,6 @@ static const char *param_var_types[PARAM_VAR_NAME_COUNT] = {
     "double", "double", "double", "double",
     "double", "double", "double", "double",
     "double","double","double", "double",
-    "double",
     "int"
 };
 //----------------------------------------------
@@ -809,10 +808,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model)
 	  continue;
         }
         if (strcmp(param_key, "Kinf_nash_surface") == 0) {
-	  //model->nash_surface_params.K_infiltration = strtod(param_value, NULL);
-	  // Kinf = C0 + C1 * S, so if Kinf is provided in the config file, set C0 and C1 as folows:
-	  model->nash_surface_params.Kinf_c0 = strtod(param_value, NULL);
-	  model->nash_surface_params.Kinf_c1 = 0.0;
+	  model->nash_surface_params.K_infiltration = strtod(param_value, NULL);
 	  is_K_infiltration_nash_surface_set = TRUE;
 	  continue;
         }
@@ -1110,10 +1106,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model)
 #if CFE_DEBUG >= 1
 	printf("Config param 'Kinf_nash_surface' not found in config file, default value is 0.05 [1/hr] \n");
 #endif
-	// model->nash_surface_params.K_infiltration  = 0.05;   // used in the runon infiltration
-	// Kinf = C0 + C1 * S, so if Kinf is not provided the default values for C0 and C1 are set below,
-	model->nash_surface_params.Kinf_c0  = 0.05;      // set C0 = value
-	model->nash_surface_params.Kinf_c1  = 0.0;	 // set C1 = 0.0
+	model->nash_surface_params.K_infiltration  = 0.05;   // used in the runon infiltration
       }
       if (is_retention_depth_nash_surface_set == FALSE) {
 #if CFE_DEBUG >= 1
@@ -1919,21 +1912,13 @@ static int Get_value_ptr (Bmi *self, const char *name, void **dest)
         return BMI_SUCCESS;
     }
 
-    if (strcmp (name, "Kinf_c0") == 0) {
+    if (strcmp (name, "K_infiltraton") == 0) {
         cfe_state_struct *cfe_ptr;
         cfe_ptr = (cfe_state_struct *) self->data;
-        *dest = (void*)&cfe_ptr->nash_surface_params.Kinf_c0;
+        *dest = (void*)&cfe_ptr->nash_surface_params.K_infiltration;
         return BMI_SUCCESS;
     }
 
-    if (strcmp (name, "Kinf_c1") == 0) {
-        cfe_state_struct *cfe_ptr;
-        cfe_ptr = (cfe_state_struct *) self->data;
-        *dest = (void*)&cfe_ptr->nash_surface_params.Kinf_c1;
-        return BMI_SUCCESS;
-    }
-      
-
     /***********************************************************/
     /***********    OUTPUT   ***********************************/
     /***********************************************************/
diff --git a/src/cfe.c b/src/cfe.c
index 100ca30..7f83551 100644
--- a/src/cfe.c
+++ b/src/cfe.c
@@ -244,10 +244,12 @@ extern void cfe(
     direct_runoff_m = giuh_convolution_integral(infiltration_excess_m,num_giuh_ordinates,
                                                 giuh_ordinates_arr,runoff_queue_m_per_timestep_arr);
   else if (surface_runoff_scheme == NASH_CASCADE)
-    direct_runoff_m =  nash_cascade_surface_runoff(infiltration_excess_m, nash_surface_params);
-
+    direct_runoff_m =  nash_cascade_surface_runoff(infiltration_excess_m, soil_reservoir_storage_deficit_m,
+						   nash_surface_params);
+  
   soil_reservoir_struct->storage_m += nash_surface_params->runon_infiltration;  // put the runon infiltrated water in the soil.
-
+  soil_reservoir_storage_deficit_m -= nash_surface_params->runon_infiltration; // update soil deficit
+  
   massbal_struct->vol_out_surface  += direct_runoff_m;
   massbal_struct->volout           += direct_runoff_m;
   massbal_struct->volout           += flux_from_deep_gw_to_chan_m;
diff --git a/src/nash_cascade.c b/src/nash_cascade.c
index b157223..d03a69d 100644
--- a/src/nash_cascade.c
+++ b/src/nash_cascade.c
@@ -40,7 +40,8 @@ extern double nash_cascade(double flux_lat_m,int num_lateral_flow_nash_reservoir
 //##############################################################
 //##############   NASH CASCADE SURFACE RUNOFF  ################
 //##############################################################
-double nash_cascade_surface_runoff(double runoff_m, struct nash_cascade_parameters *nash_params)
+double nash_cascade_surface_runoff(double runoff_m, double soil_storage_deficit_m,
+				   struct nash_cascade_parameters *nash_params)
 {
   //##############################################################
   // Solve for the flow through the Nash cascade to delay the
@@ -51,9 +52,7 @@ double nash_cascade_surface_runoff(double runoff_m, struct nash_cascade_paramete
   int N_nash    = nash_params->N_nash;
   double K_nash = nash_params->K_nash;
   double depth_r = nash_params->retention_depth;
-  //double K_infil = nash_params->K_infiltration;
-  double c0 = nash_params->Kinf_c0;
-  double c1 = nash_params->Kinf_c1;
+  double K_infil = nash_params->K_infiltration;
   
   // local vars
   double dt_h  = 1.0;             // model timestep [hour]
@@ -64,10 +63,11 @@ double nash_cascade_surface_runoff(double runoff_m, struct nash_cascade_paramete
   double Q_r;                    // discharge from reservoir
   double Q_out;                  // discharge at the outlet (the last reservoir) per subtimestep
   double Q_infil;                // discharge from reservoirs to soil
-  double K_infil;                  // K infiltration = c0 + c1 * S
   double Q_to_channel_m = 0.0;  // total outflow to channel per timestep
   double Q_to_soil_m    = 0.0;  // runon infiltration (losses from surface runoff to soil)
-
+  double soil_deficit_m = soil_storage_deficit_m; // local variable to track the soil storage deficit
+  double infil_m;
+  
   nash_params->nash_storage[0] += runoff_m;
 
   // Loop through number of sub-timesteps
@@ -76,9 +76,44 @@ double nash_cascade_surface_runoff(double runoff_m, struct nash_cascade_paramete
     //Loop through reservoirs
     for(int i = 0; i < N_nash; i++) {
 
+      // First: infiltration capacity should be satisfied before routing water through Nash reservoirs
+
+      // compute runon infiltration
+      Q_infil = K_infil * nash_params->nash_storage[i];
+      infil_m = Q_infil * subdt;
+
+      // case 1: soil_deficit greater than infil_m, all water can infiltrate, and soil_deficit decreases
+      // case 2: soil_deficit less than infil_m, portion of water can infiltrate, soil_deficit gets zero
+      // case 3: soil_deficit is zero, no infiltration
+
+      // determine the right amount of infiltrated water based on soil deficit
+      if (soil_deficit_m > 0.0 && soil_deficit_m > infil_m) {
+	soil_deficit_m -= infil_m;           // update soil deficit for the next subtimstep iteration
+      }
+      else if (soil_deficit_m > 0.0) {
+	infil_m = soil_deficit_m;
+	soil_deficit_m = 0.0; 
+      }
+      else {
+	infil_m = 0.0; // if got here, then soil is fully saturated, so no infiltration
+      }
+
+      // update nash storage (subtract infiltrated water)
+      if  (nash_params->nash_storage[i] >= infil_m)
+        nash_params->nash_storage[i] = nash_params->nash_storage[i] - infil_m;
+      else {
+        infil_m = nash_params->nash_storage[i];
+        nash_params->nash_storage[i] = 0.0;
+      }
+
+
+      Q_to_soil_m += infil_m;  // water volume that infiltrates per subtimestep from each reservoir
+
+      /*=========================================================================================*/
+      // Second: time to route water through Nash reservoirs
       S = nash_params->nash_storage[i];
 
-      // first reservoir
+      // determine the amount of surface water available for routing from the first reservoir
       if (i == 0 &&  S > depth_r)
         S -= depth_r;
       else if (i == 0)
@@ -92,21 +127,7 @@ double nash_cascade_surface_runoff(double runoff_m, struct nash_cascade_paramete
         nash_params->nash_storage[i+1] += dS;
       else
         Q_out = Q_r;
-
-      //compute runon infiltration
-      K_infil = c0 + c1 * nash_params->nash_storage[i];
-      Q_infil = K_infil * nash_params->nash_storage[i];
-      dS_infil = Q_infil * subdt;
-
-      if  (nash_params->nash_storage[i] >= dS_infil)
-        nash_params->nash_storage[i] = nash_params->nash_storage[i] - dS_infil;
-      else {
-        dS_infil = nash_params->nash_storage[i];
-        nash_params->nash_storage[i] = 0.0;
-      }
-
-      Q_to_soil_m += dS_infil;  // water volume that infiltrates per subtimestep from each reservoir
-
+      
     }
 
    Q_to_channel_m += Q_out * subdt;  // Q_r at the end of N_nash loop is the discharge at the outlet