diff --git a/Phreeqc.cpp b/Phreeqc.cpp
index d89c6dec2..25d903afa 100644
--- a/Phreeqc.cpp
+++ b/Phreeqc.cpp
@@ -1722,7 +1722,10 @@ Phreeqc::InternalCopy(const Phreeqc* pSrc)
 	viscos_0 = pSrc->viscos_0;
 	viscos_0_25 = pSrc->viscos_0_25; // viscosity of the solution, of pure water, of pure water at 25 C
 	density_x = pSrc->density_x;
-	cell_pore_volume = pSrc->cell_pore_volume;;
+	solution_volume_x = pSrc->solution_volume_x;
+	solution_mass_x = pSrc->solution_mass_x;
+	kgw_kgs = pSrc->kgw_kgs;
+	cell_pore_volume = pSrc->cell_pore_volume;
 	cell_porosity = pSrc->cell_porosity;
 	cell_volume = pSrc->cell_volume;
 	cell_saturation = pSrc->cell_saturation;
diff --git a/Phreeqc.h b/Phreeqc.h
index 225c52bd5..be150845f 100644
--- a/Phreeqc.h
+++ b/Phreeqc.h
@@ -1523,6 +1523,7 @@ class Phreeqc
 	int iterations;
 	int gamma_iterations;
 	size_t density_iterations;
+	LDBLE kgw_kgs;
 	int run_reactions_iterations;
 	int overall_iterations;
 
@@ -1629,6 +1630,8 @@ class Phreeqc
 	int print_viscosity;
 	LDBLE viscos, viscos_0, viscos_0_25; // viscosity of the solution, of pure water, of pure water at 25 C
 	LDBLE density_x;
+	LDBLE solution_volume_x;
+	LDBLE solution_mass_x;
 	LDBLE cell_pore_volume;
 	LDBLE cell_porosity;
 	LDBLE cell_volume;
diff --git a/basicsubs.cpp b/basicsubs.cpp
index 7b5a6f285..14744601a 100644
--- a/basicsubs.cpp
+++ b/basicsubs.cpp
@@ -592,6 +592,8 @@ calc_dens(void)
 	{
 		density_x = rho_0 * (1e3 + M_T / mass_water_aq_x) / (rho_0 * V_solutes / mass_water_aq_x + 1e3);
 	}
+	solution_mass_x = 1e-3*(M_T + s_h2o->moles * s_h2o->gfw);
+	solution_volume_x = solution_mass_x / density_x;
 	return density_x;
 	//M_T /= 1e3;
 	//solution_mass =  mass_water_aq_x + M_T;
@@ -604,7 +606,7 @@ calc_dens(void)
 }
 /* VP: Density End */
 /* DP: Function for interval halving */
-
+#ifdef NOT_USED
 LDBLE Phreeqc::
 f_rho(LDBLE rho_old, void* cookie)
 /* ---------------------------------------------------------------------- */
@@ -622,7 +624,8 @@ f_rho(LDBLE rho_old, void* cookie)
 	rho = rho + pThis->rho_0;
 	return (rho - rho_old);
 }
-
+#endif
+#ifdef original
 /* ---------------------------------------------------------------------- */
 LDBLE Phreeqc::
 calc_solution_volume(void)
@@ -653,7 +656,19 @@ calc_solution_volume(void)
 	LDBLE vol = 1e-3 * total_mass / rho;
 	return (vol);
 }
-
+#endif
+/* ---------------------------------------------------------------------- */
+LDBLE Phreeqc::
+calc_solution_volume(void)
+/* ---------------------------------------------------------------------- */
+{
+	/*
+	 *   Calculates solution volume based on sum of mass of element plus density
+	 */
+	LDBLE rho = calc_dens();
+	LDBLE vol = solution_volume_x;
+	return (vol);
+}
 /* ---------------------------------------------------------------------- */
 LDBLE Phreeqc::
 calc_logk_n(const char* name)
diff --git a/mainsubs.cpp b/mainsubs.cpp
index e57c6389c..051e85be0 100644
--- a/mainsubs.cpp
+++ b/mainsubs.cpp
@@ -409,6 +409,8 @@ initial_solutions(int print)
 					set(TRUE);
 					converge = model();
 				}
+				calc_dens();
+				kgw_kgs = mass_water_aq_x / solution_mass_x;
 				density_iterations++;
 				if (solution_ref.Get_initial_data()->Get_calc_density())
 				{
diff --git a/prep.cpp b/prep.cpp
index e08cd5f8f..7494dfbed 100644
--- a/prep.cpp
+++ b/prep.cpp
@@ -1942,6 +1942,10 @@ convert_units(cxxSolution *solution_ptr)
 		strstr(initial_data_ptr->Get_units().c_str(), "/l") != NULL)
 	{
 		mass_water_aq_x = 1.0 - 1e-3 * sum_solutes;
+		if (density_iterations > 0)
+		{
+			mass_water_aq_x = kgw_kgs;
+		}
 		if (mass_water_aq_x <= 0)
 		{
 			error_string = sformatf( "Solute mass exceeds solution mass in conversion from /kgs to /kgw.\n"