diff --git a/phreeqcpp/basicsubs.cpp b/phreeqcpp/basicsubs.cpp
index 003726aa..c48d56b7 100644
--- a/phreeqcpp/basicsubs.cpp
+++ b/phreeqcpp/basicsubs.cpp
@@ -2445,46 +2445,61 @@ surf_total(const char *total_name, const char *surface_name)
 
 		// surface matches, now match element or redox state
 		struct rxn_token *rxn_ptr;
-		for (rxn_ptr = s_x[j]->rxn_s->token + 1; rxn_ptr->s != NULL; rxn_ptr++)
+		if (s_x[j]->mole_balance == NULL)
 		{
-			if (redox && rxn_ptr->s->secondary)
+			for (rxn_ptr = s_x[j]->rxn_s->token + 1; rxn_ptr->s != NULL; rxn_ptr++)
 			{
-				token = rxn_ptr->s->secondary->elt->name;
-			}
-			else if (!redox && rxn_ptr->s->secondary)
-			{
-				token = rxn_ptr->s->secondary->elt->primary->elt->name;
-			}
-			else if (!redox && rxn_ptr->s->primary)
-			{
-				token = rxn_ptr->s->primary->elt->name;
-			}
-			else
-			{
-				continue;
-			}
-			if (strcmp(token.c_str(), total_name) == 0)
-			{
-				t += rxn_ptr->coef * s_x[j]->moles;
-				break;
-			}
-			else
-			// sum all sites in case total_name is a surface name without underscore surf ("Hfo_w", "Hfo")
-			{
-				if (rxn_ptr->s->type == SURF)
+				if (redox && rxn_ptr->s->secondary)
 				{
-					if (token.find("_") != std::string::npos)
+					token = rxn_ptr->s->secondary->elt->name;
+				}
+				else if (!redox && rxn_ptr->s->secondary)
+				{
+					token = rxn_ptr->s->secondary->elt->primary->elt->name;
+				}
+				else if (!redox && rxn_ptr->s->primary)
+				{
+					token = rxn_ptr->s->primary->elt->name;
+				}
+				else
+				{
+					continue;
+				}
+				if (strcmp(token.c_str(), total_name) == 0)
+				{
+					t += rxn_ptr->coef * s_x[j]->moles;
+					break;
+				}
+				else
+					// sum all sites in case total_name is a surface name without underscore surf ("Hfo_w", "Hfo")
+				{
+					if (rxn_ptr->s->type == SURF)
 					{
-						token = token.substr(0, token.find("_"));
-					}
-					if (strcmp(token.c_str(), total_name) == 0)
-					{
-						t += rxn_ptr->coef * s_x[j]->moles;
-						break;
+						if (token.find("_") != std::string::npos)
+						{
+							token = token.substr(0, token.find("_"));
+						}
+						if (strcmp(token.c_str(), total_name) == 0)
+						{
+							t += rxn_ptr->coef * s_x[j]->moles;
+							break;
+						}
 					}
 				}
 			}
 		}
+		else
+		{
+			for (int i = 0; s_x[j]->next_secondary[i].elt != NULL; i++)
+			{
+				token = s_x[j]->next_secondary[i].elt->name;
+				if (strcmp(token.c_str(), total_name) == 0)
+				{
+					t += s_x[j]->next_secondary[i].coef * s_x[j]->moles;
+					break;
+				}
+			}
+		}
 	}
 	return t;
 }
diff --git a/phreeqcpp/model.cpp b/phreeqcpp/model.cpp
index f89990a2..3086158c 100644
--- a/phreeqcpp/model.cpp
+++ b/phreeqcpp/model.cpp
@@ -532,7 +532,7 @@ check_residuals(void)
 		{
 			if (x[i]->ss_in == FALSE)
 				continue;
-			if (x[i]->moles <= MIN_TOTAL_SS)
+			if (x[i]->moles <= 1e2*MIN_TOTAL)
 				continue;
 			if (residual[i] >= epsilon
 				|| residual[i] <= -epsilon /* || stop_program == TRUE */ )
@@ -1597,7 +1597,9 @@ ineq(int in_kode)
 					   (size_t) (count_unknowns + 1) * sizeof(LDBLE));
 				ineq_array[l_count_rows * max_column_count + i] = 1.0;
 				ineq_array[l_count_rows * max_column_count + count_unknowns] =
-					0.99 * x[i]->moles - MIN_TOTAL_SS;
+					x[i]->moles;
+					//0.99 * x[i]->moles - MIN_TOTAL_SS;
+					//0.9 * x[i]->moles;
 				back_eq[l_count_rows] = i;
 				l_count_rows++;
 			}
@@ -3264,6 +3266,20 @@ reset(void)
 						}
 						factor = f0;
 					}
+					if (x[i]->moles > 0 && -delta[i] > 1e1*x[i]->moles)
+					{
+						f0 = (-delta[i])/(1e1*x[i]->moles);
+						if (f0 > factor)
+						{
+							if (debug_model == TRUE)
+							{
+								output_msg(sformatf(
+									"%-10.10s, Precipitating too much mineral.\t%f\n",
+									x[i]->description, (double)f0));
+							}
+							factor = f0;
+						}
+					}
 				}
 			}
 		}
@@ -4207,6 +4223,7 @@ residuals(void)
 		else if (x[i]->type == SS_MOLES)
 		{
 			residual[i] = x[i]->f * LOG_10;
+			if (x[i]->moles <= 1e2*MIN_TOTAL) continue;
 			if (fabs(residual[i]) > l_toler && x[i]->ss_in == TRUE)
 			{
 				if (print_fail)
diff --git a/phreeqcpp/transport.cpp b/phreeqcpp/transport.cpp
index d4eb8df1..0f573de4 100644
--- a/phreeqcpp/transport.cpp
+++ b/phreeqcpp/transport.cpp
@@ -7,6 +7,7 @@
 #include "SSassemblage.h"
 #include "cxxKinetics.h"
 #include "Solution.h"
+#include <limits>
 
 LDBLE F_Re3 = F_C_MOL / (R_KJ_DEG_MOL * 1e3);
 LDBLE tk_x2; // average tx_x of icell and jcell
@@ -1107,8 +1108,15 @@ init_mix(void)
 		}
 		else
 		{
+			const int max_int = (int) (std::numeric_limits<int>::max)();
+			if ((int)round(2.25 * maxmix) > max_int)
+			{
+				char token[MAX_LENGTH];
+				sprintf(token, "Calculated number of mixes %g, is beyond program limit,\nERROR: please decrease time_step.", 2.25 * maxmix);
+				error_msg(token, STOP);
+			}
 			if (bcon_first == 1 || bcon_last == 1)
-				l_nmix = 1 + (int) floor(2.25 * maxmix);
+				l_nmix = 1 + (long int) floorl(2.25 * maxmix);
 			else
 				l_nmix = 1 + (int) floor(1.5 * maxmix);
 
@@ -2046,6 +2054,7 @@ multi_D(LDBLE DDt, int mobile_cell, int stagnant)
 				* 1. obtain J_ij...
 				*/
 				il_calcs = find_J(icell, jcell, mixf, DDt, stagnant);
+
 				if (find_current)
 				{
 					if (i < last_c)
@@ -2078,6 +2087,8 @@ multi_D(LDBLE DDt, int mobile_cell, int stagnant)
 						continue;
 					}
 				}
+				if (!ct[icell].J_ij_count_spec)
+					continue;
 
 				/*
 				* 2. sum up the primary or secondary master_species