diff --git a/CHANGELOG.md b/CHANGELOG.md
index e3f64abb..c5915689 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,6 +13,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 
+- `Solution`: Fixed an issue in which repeated calls to `equilibrate` when using `NativeEOS` or `PHREEQCEOS` would
+  change the mass of the `Solution` slightly. This was attributed to the fact that `pyEQL` and `PHREEQC` use slightly
+  different molecular weights for water.
 - `Solution`: `get_total_amount` and related methods could fail when the oxidation state of an element was
   unknown (e.g., 'Br') (Issue [#116](https://github.com/KingsburyLab/pyEQL/issues/116))
 
diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py
index bb645f58..30a30b5c 100644
--- a/src/pyEQL/engines.py
+++ b/src/pyEQL/engines.py
@@ -5,6 +5,7 @@
 :license: LGPL, see LICENSE for more details.
 
 """
+
 import os
 import warnings
 from abc import ABC, abstractmethod
@@ -189,7 +190,6 @@ def _setup_ppsol(self, solution):
             # PHREEQC will assume 1 kg if not specified, there is also no direct way to specify volume, so we
             # really have to specify the solvent mass in 1 liter of solution
             "water": solv_mass,
-            "density": solution.density.to("g/mL").magnitude,
         }
         if solution.balance_charge == "pH":
             d["pH"] = str(d["pH"]) + " charge"
@@ -234,9 +234,7 @@ def _setup_ppsol(self, solution):
         self.ppsol = ppsol
 
     def _destroy_ppsol(self):
-        """
-        Remove the PhreeqPython solution from memory
-        """
+        """Remove the PhreeqPython solution from memory"""
         if self.ppsol is not None:
             self.ppsol.forget()
             self.ppsol = None
@@ -370,9 +368,7 @@ def get_activity_coefficient(self, solution, solute):
             )
 
             logger.info(
-                "Calculated activity coefficient of species {} as {} based on salt {} using Pitzer model".format(
-                    solute, activity_coefficient, salt
-                )
+                f"Calculated activity coefficient of species {solute} as {activity_coefficient} based on salt {salt} using Pitzer model"
             )
             molal = activity_coefficient
 
@@ -549,9 +545,7 @@ def get_osmotic_coefficient(self, solution):
                 )
 
                 logger.info(
-                    "Calculated osmotic coefficient of water as {} based on salt {} using Pitzer model".format(
-                        osmotic_coefficient, item.formula
-                    )
+                    f"Calculated osmotic coefficient of water as {osmotic_coefficient} based on salt {item.formula} using Pitzer model"
                 )
                 effective_osmotic_sum += concentration * osmotic_coefficient
 
@@ -658,7 +652,9 @@ def equilibrate(self, solution):
         if self.ppsol is not None:
             self.ppsol.forget()
         self._setup_ppsol(solution)
-        # self._stored_comp = solution.components
+
+        # store the original solvent mass
+        orig_solvent_moles = solution.components[solution.solvent]
 
         # use the output from PHREEQC to update the Solution composition
         # the .species_moles attribute should return MOLES (not moles per ___)
@@ -679,17 +675,28 @@ def equilibrate(self, solution):
             )
         for s in missing_species:
             charge_adjust += -1 * solution.get_amount(s, "eq").magnitude
+        if charge_adjust != 0:
+            logger.warning(
+                "After equilibration, the charge balance of the solution was not electroneutral."
+                f" {charge_adjust} eq of charge were added via {solution.balance_charge}"
+            )
 
-        # re-adjust charge balance
-        if solution.balance_charge is None:
-            pass
-        elif solution.balance_charge == "pH":
-            solution.components["H+"] += charge_adjust
-        elif solution.balance_charge == "pE":
-            raise NotImplementedError
-        else:
-            z = solution.get_property(solution.balance_charge, "charge")
-            solution.add_amount(solution.balance_charge, f"{charge_adjust/z} mol")
+            # re-adjust charge balance
+            if solution.balance_charge is None:
+                pass
+            elif solution.balance_charge == "pH":
+                solution.components["H+"] += charge_adjust
+            elif solution.balance_charge == "pE":
+                raise NotImplementedError
+            else:
+                z = solution.get_property(solution.balance_charge, "charge")
+                solution.add_amount(solution.balance_charge, f"{charge_adjust/z} mol")
+
+        # rescale the solvent mass to ensure the total mass of solution does not change
+        # this is important because PHREEQC and the pyEQL database may use slightly different molecular
+        # weights for water. Since water amount is passed to PHREEQC in kg but returned in moles, each
+        # call to equilibrate can thus result in a slight change in the Solution mass.
+        solution.components[solution.solvent] = orig_solvent_moles
 
     def __deepcopy__(self, memo):
         # custom deepcopy required because the PhreeqPython instance used by the Native and Phreeqc engines
diff --git a/tests/test_phreeqc.py b/tests/test_phreeqc.py
index 20850dce..835b000f 100644
--- a/tests/test_phreeqc.py
+++ b/tests/test_phreeqc.py
@@ -168,18 +168,25 @@ def test_equilibrate(s1, s2, s5_pH, s6_Ca, caplog):
     assert "H2(aq)" not in s1.components
     orig_pH = s1.pH
     orig_pE = s1.pE
+    orig_mass = s1.mass
+    orig_density = s2.density.magnitude
+    orig_solv_mass = s2.solvent_mass.magnitude
     s1.equilibrate()
     assert "H2(aq)" in s1.components
     assert np.isclose(s1.charge_balance, 0, atol=1e-8)
     assert np.isclose(s1.pH, orig_pH, atol=0.01)
     assert np.isclose(s1.pE, orig_pE)
+    # assert np.isclose(s1.mass, orig_mass)
+    # assert np.isclose(s1.density.magnitude, orig_density)
+    # assert np.isclose(s1.solvent_mass.magnitude, orig_solv_mass)
 
     assert "NaOH(aq)" not in s2.components
-    s2.equilibrate()
-    orig_pH = s2.pH
-    orig_pE = s2.pE
     orig_density = s2.density.magnitude
     orig_solv_mass = s2.solvent_mass.magnitude
+    orig_pH = s2.pH
+    orig_pE = s2.pE
+    orig_mass = s2.mass
+    s2.equilibrate()
     assert "NaOH(aq)" in s2.components
 
     # total element concentrations should be conserved after equilibrating
@@ -187,14 +194,18 @@ def test_equilibrate(s1, s2, s5_pH, s6_Ca, caplog):
     assert np.isclose(s2.get_total_amount("Cl", "mol").magnitude, 8)
     assert np.isclose(s2.solvent_mass.magnitude, orig_solv_mass)
     assert np.isclose(s2.density.magnitude, orig_density)
+    # the pH should drop due to hydrolysis of Ca+2
+    assert s2.pH < orig_pH
+    assert np.isclose(s2.pE, orig_pE)
+    assert np.isclose(s2.mass, orig_mass)
     # this solution has balance_charge=None, therefore, the charge balance
     # may be off after equilibration
     assert not np.isclose(s2.charge_balance, 0, atol=1e-8)
-    assert np.isclose(s2.pH, orig_pH, atol=0.01)
-    assert np.isclose(s2.pE, orig_pE)
+    eq_Hplus = s2.components["H+"]
     s2.balance_charge = "pH"
     s2.equilibrate()
     assert np.isclose(s2.charge_balance, 0, atol=1e-8)
+    assert s2.components["H+"] > eq_Hplus
 
     # test log message if there is a species not present in the phreeqc database
     s_zr = Solution({"Zr+4": "0.05 mol/kg", "Na+": "0.05 mol/kg", "Cl-": "0.1 mol/kg"}, engine="phreeqc")