diff --git a/requirements-dev.txt b/requirements-dev.txt
index 6f633fe494..8033542896 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -13,4 +13,5 @@ sphinx==7.1.*
 nbsphinx
 sphinx_rtd_theme
 
+idaes-pse @ git+https://github.com/IDAES/idaes-pse.git@2.10.0rc0
 -e .[testing,notebooks,oli_api]
diff --git a/setup.py b/setup.py
index e9a1233d37..77a894be5a 100644
--- a/setup.py
+++ b/setup.py
@@ -69,7 +69,7 @@
         # primary requirements for unit and property models
         # allow X.Y.Z stable release(s) and X.Y+1.dev0 (i.e. the main branch after X.Y.Z)
         # disallow X.Y+1.0rc0 (i.e. forcing a manual update to this requirement)
-        "idaes-pse >=2.9.0,<2.10.0rc0",
+        "idaes-pse @ git+https://github.com/IDAES/idaes-pse.git@2.10.0rc0",
         "pyomo>=6.6.1",
         "watertap-solvers",
         "pyyaml",  # watertap.core.wt_database
diff --git a/watertap/core/membrane_channel1d.py b/watertap/core/membrane_channel1d.py
index a0269173fa..0e998e1d96 100644
--- a/watertap/core/membrane_channel1d.py
+++ b/watertap/core/membrane_channel1d.py
@@ -242,14 +242,14 @@ def initialize(
             self.release_state(source_flags, outlvl)
 
     def calculate_scaling_factors(self):
+        if iscale.get_scaling_factor(self.area) is None:
+            iscale.set_scaling_factor(self.area, 100)
+
         super().calculate_scaling_factors()
 
         # setting scaling factors for variables
 
         # will not override if the user provides the scaling factor
-        ## default of 1 set by ControlVolume1D
-        if iscale.get_scaling_factor(self.area) == 1:
-            iscale.set_scaling_factor(self.area, 100)
 
         if hasattr(self, "pressure_change_total"):
             for v in self.pressure_change_total.values():
diff --git a/watertap/core/membrane_channel_base.py b/watertap/core/membrane_channel_base.py
index 0b3f365599..fb4fd89a1c 100644
--- a/watertap/core/membrane_channel_base.py
+++ b/watertap/core/membrane_channel_base.py
@@ -19,6 +19,7 @@
     Expression,
     Param,
     NonNegativeReals,
+    value,
     Var,
     units as pyunits,
 )
@@ -996,7 +997,11 @@ def calculate_scaling_factors(self):
         if hasattr(self, "velocity"):
             for v in self.velocity.values():
                 if iscale.get_scaling_factor(v) is None:
-                    iscale.set_scaling_factor(v, 1)
+                    sf_A = iscale.get_scaling_factor(self.area, default=1)
+                    sf_F_vol = iscale.get_scaling_factor(
+                        self.properties[t, x].flow_vol_phase["Liq"], default=1
+                    )
+                    iscale.set_scaling_factor(v, sf_F_vol / sf_A, 1)
 
         if hasattr(self, "friction_factor_darcy"):
             for v in self.friction_factor_darcy.values():
@@ -1008,6 +1013,71 @@ def calculate_scaling_factors(self):
                 if iscale.get_scaling_factor(v) is None:
                     iscale.set_scaling_factor(v, 1)
 
+        # TODO why did this file have zero scaling factors for constraints
+        if hasattr(self, "eq_dh"):
+            sf_dh = iscale.get_scaling_factor(self.dh, default=1)
+            iscale.constraint_scaling_transform(self.eq_dh, sf_dh)
+
+        if hasattr(self, "eq_area"):
+            sf_A = iscale.get_scaling_factor(self.area, default=1)
+            iscale.constraint_scaling_transform(self.eq_area, sf_A)
+
+        if hasattr(self, "eq_velocity"):
+            for (t, x), condata in self.eq_velocity.items():
+                sf_F_vol = iscale.get_scaling_factor(
+                    self.properties[t, x].flow_vol_phase["Liq"], default=1
+                )
+                iscale.constraint_scaling_transform(condata, sf_F_vol)
+
+        if hasattr(self, "eq_equal_flow_vol_interface"):
+            for (t, x), condata in self.eq_equal_flow_vol_interface.items():
+                sf_F_vol = min(
+                    iscale.get_scaling_factor(
+                        self.properties[t, x].flow_vol_phase["Liq"], default=1
+                    ),
+                    iscale.get_scaling_factor(
+                        self.properties_interface[t, x].flow_vol_phase["Liq"], default=1
+                    ),
+                )
+                iscale.constraint_scaling_transform(condata, sf_F_vol)
+
+        if hasattr(self, "eq_equal_pressure_interface"):
+            for (t, x), condata in self.eq_equal_pressure_interface.items():
+                sf_P = min(
+                    iscale.get_scaling_factor(
+                        self.properties_interface[t, x].pressure, default=1
+                    ),
+                    iscale.get_scaling_factor(
+                        self.properties[t, x].pressure, default=1
+                    ),
+                )
+                iscale.constraint_scaling_transform(condata, sf_P)
+
+        if hasattr(self, "eq_K"):
+            for (t, x, j), condata in self.eq_K.items():
+                sf_K = iscale.get_scaling_factor(self.K[t, x, j])
+                # Hopefully h has been calculated before this
+                sf_h = 1 / value(self.dh)
+                iscale.constraint_scaling_transform(condata, sf_K * sf_h)
+
+        if hasattr(self, "eq_N_Re"):
+            for (t, x), condata in self.eq_N_Re.items():
+                sf_Re = iscale.get_scaling_factor(self.N_Re[t, x], default=1)
+                sf_A = iscale.get_scaling_factor(
+                    self.area, default=1 / value(self.area)
+                )
+                sf_visc = iscale.get_scaling_factor(
+                    self.properties[t, x].visc_d_phase["Liq"], default=1
+                )
+                iscale.constraint_scaling_transform(condata, sf_Re * sf_A * sf_visc)
+
+        if hasattr(self, "eq_N_Sc_comp"):
+            for (t, x, j), condata in self.eq_N_Sc_comp.items():
+                sf_visc = iscale.get_scaling_factor(
+                    self.properties[t, x].visc_d_phase["Liq"], default=1
+                )
+                iscale.constraint_scaling_transform(condata, sf_visc)
+
 
 # helper for validating configuration arguments for this CV
 def validate_membrane_config_args(unit):
diff --git a/watertap/core/zero_order_sido.py b/watertap/core/zero_order_sido.py
index 3f0ecf1bcb..43ce41f21c 100644
--- a/watertap/core/zero_order_sido.py
+++ b/watertap/core/zero_order_sido.py
@@ -282,16 +282,17 @@ def initialize_sido(
 
 def calculate_scaling_factors_sido(self):
     # Get default scale factors and do calculations from base classes
-    for t, v in self.water_recovery_equation.items():
-        iscale.constraint_scaling_transform(
-            v,
-            iscale.get_scaling_factor(
-                self.properties_in[t].flow_mass_comp["H2O"],
-                default=1,
-                warning=True,
-                hint=" for water recovery",
-            ),
-        )
+    if hasattr(self, "water_recovery_equation"):
+        for t, v in self.water_recovery_equation.items():
+            iscale.constraint_scaling_transform(
+                v,
+                iscale.get_scaling_factor(
+                    self.properties_in[t].flow_mass_comp["H2O"],
+                    default=1,
+                    warning=True,
+                    hint=" for water recovery",
+                ),
+            )
 
     for t, v in self.water_balance.items():
         iscale.constraint_scaling_transform(
@@ -320,6 +321,15 @@ def calculate_scaling_factors_sido(self):
             ),
         )  # would just be a duplicate of above
 
+    if hasattr(self, "eq_isobaric"):
+        for (t, port), condata in self.eq_isobaric.items():
+            obj = getattr(self, port)
+            sf_P = min(
+                iscale.get_scaling_factor(self.inlet.pressure[t], default=1),
+                iscale.get_scaling_factor(obj.pressure[t], default=1),
+            )
+            iscale.constraint_scaling_transform(condata, sf_P)
+
 
 def _get_Q_sido(self, t):
     return self.properties_in[t].flow_vol
diff --git a/watertap/flowsheets/full_water_resource_recovery_facility/BSM2.py b/watertap/flowsheets/full_water_resource_recovery_facility/BSM2.py
index 623c8a6b40..38c21fe160 100644
--- a/watertap/flowsheets/full_water_resource_recovery_facility/BSM2.py
+++ b/watertap/flowsheets/full_water_resource_recovery_facility/BSM2.py
@@ -1042,7 +1042,9 @@ def initialize_system(m):
     seq.set_guesses_for(m.fs.R1.inlet, tear_guesses1)
     seq.set_guesses_for(m.fs.asm_adm.inlet, tear_guesses2)
 
-    initializer = BlockTriangularizationInitializer()
+    initializer = BlockTriangularizationInitializer(
+        calculate_variable_options={"eps": 2e-8}, skip_final_solve=True
+    )
 
     def function(unit):
         try:
diff --git a/watertap/flowsheets/full_water_resource_recovery_facility/BSM2_P_extension.py b/watertap/flowsheets/full_water_resource_recovery_facility/BSM2_P_extension.py
index 66e83bb0e9..72f6611345 100644
--- a/watertap/flowsheets/full_water_resource_recovery_facility/BSM2_P_extension.py
+++ b/watertap/flowsheets/full_water_resource_recovery_facility/BSM2_P_extension.py
@@ -502,14 +502,20 @@ def scale_system(m, bio_P=False):
     ad_scaler = ADScaler()
     ad_scaler.scale_model(m.fs.AD)
 
-    set_scaling_factor(m.fs.AD.KH_co2, 1e1)
-    set_scaling_factor(m.fs.AD.KH_ch4, 1e1)
-    set_scaling_factor(m.fs.AD.KH_h2, 1e2)
-    set_scaling_factor(m.fs.AD.liquid_phase.heat, 1e1)
+    for vardata in m.fs.AD.KH_co2.values():
+        set_scaling_factor(vardata, 1e1)
+    for vardata in m.fs.AD.KH_ch4.values():
+        set_scaling_factor(vardata, 1e1)
+    for vardata in m.fs.AD.KH_h2.values():
+        set_scaling_factor(vardata, 1e2)
+    for vardata in m.fs.AD.liquid_phase.heat.values():
+        set_scaling_factor(vardata, 1e1)
     if bio_P:
-        set_scaling_factor(m.fs.AD.liquid_phase.reactions[0].S_H, 1e1)
+        for blkdata in m.fs.AD.liquid_phase.reactions.values():
+            set_scaling_factor(blkdata.S_H, 1e1)
     else:
-        set_scaling_factor(m.fs.AD.liquid_phase.reactions[0].S_H, 1e2)
+        for blkdata in m.fs.AD.liquid_phase.reactions.values():
+            set_scaling_factor(blkdata.S_H, 1e2)
 
     cstr_list = [m.fs.R1, m.fs.R2, m.fs.R3, m.fs.R4]
     cstr_scaler = CSTRScaler()
@@ -517,7 +523,8 @@ def scale_system(m, bio_P=False):
         cstr_scaler.scale_model(unit)
 
     for unit in cstr_list:
-        set_scaling_factor(unit.hydraulic_retention_time, 1e-3)
+        for vardata in unit.hydraulic_retention_time.values():
+            set_scaling_factor(vardata, 1e-3)
 
     aeration_list = [m.fs.R5, m.fs.R6, m.fs.R7]
     aeration_scaler = AerationTankScaler()
@@ -525,29 +532,39 @@ def scale_system(m, bio_P=False):
         aeration_scaler.scale_model(unit)
 
     for R in aeration_list:
-        set_scaling_factor(R.KLa, 1e-1)
+        for vardata in R.KLa.values():
+            set_scaling_factor(vardata, 1e-1)
         if bio_P:
-            set_scaling_factor(R.hydraulic_retention_time[0], 1e-2)
+            for vardata in R.hydraulic_retention_time.values():
+                set_scaling_factor(vardata, 1e-2)
 
     reactor_list = [m.fs.R1, m.fs.R2, m.fs.R3, m.fs.R4, m.fs.R5, m.fs.R6, m.fs.R7]
     for unit in reactor_list:
-        set_scaling_factor(unit.control_volume.reactions[0.0].rate_expression, 1e3)
-        set_scaling_factor(unit.cstr_performance_eqn, 1e3)
-        set_scaling_factor(
-            unit.control_volume.rate_reaction_stoichiometry_constraint, 1e3
-        )
-        set_scaling_factor(unit.control_volume.material_balances, 1e3)
+        for blkdata in unit.control_volume.reactions.values():
+            for vardata in blkdata.rate_expression.values():
+                set_scaling_factor(vardata, 1e3)
+        for condata in unit.cstr_performance_eqn.values():
+            set_scaling_factor(condata, 1e3)
+        for (
+            condata
+        ) in unit.control_volume.rate_reaction_stoichiometry_constraint.values():
+            set_scaling_factor(condata, 1e3)
+        for condata in unit.control_volume.material_balances.values():
+            set_scaling_factor(condata, 1e3)
 
     if bio_P:
-        set_scaling_factor(
-            m.fs.R5.control_volume.rate_reaction_generation[0, "Liq", "S_I"], 1e-3
-        )
-        constraint_scaling_transform(
-            m.fs.R5.control_volume.rate_reaction_stoichiometry_constraint[
-                0, "Liq", "H2O"
-            ],
-            1e-6,
-        )
+        for t in m.fs.time:
+            set_scaling_factor(
+                m.fs.R5.control_volume.rate_reaction_generation[t, "Liq", "S_I"], 1e-3
+            )
+            # TODO is it intentional that we're setting scaling factors for some constraints
+            # but doing constraint scaling transforms for others?
+            constraint_scaling_transform(
+                m.fs.R5.control_volume.rate_reaction_stoichiometry_constraint[
+                    t, "Liq", "H2O"
+                ],
+                1e-6,
+            )
 
     clarifier_list = [m.fs.CL, m.fs.CL2]
     clarifier_scaler = ClarifierScaler()
@@ -566,7 +583,8 @@ def scale_system(m, bio_P=False):
     ad_as_scaler = ADM1ASM2dScaler()
     ad_as_scaler.scale_model(m.fs.translator_adm1_asm2d)
 
-    set_scaling_factor(m.fs.P1.control_volume.work[0], 1e-2)
+    for vardata in m.fs.P1.control_volume.work.values():
+        set_scaling_factor(vardata, 1e-2)
 
     for var in m.fs.component_data_objects(pyo.Var, descend_into=True):
         if "flow_vol" in var.name:
diff --git a/watertap/flowsheets/full_water_resource_recovery_facility/test/test_BSM2_P_extension.py b/watertap/flowsheets/full_water_resource_recovery_facility/test/test_BSM2_P_extension.py
index 6facc1771f..db73f37ee3 100644
--- a/watertap/flowsheets/full_water_resource_recovery_facility/test/test_BSM2_P_extension.py
+++ b/watertap/flowsheets/full_water_resource_recovery_facility/test/test_BSM2_P_extension.py
@@ -75,17 +75,9 @@ def test_structural_issues(self, system_frame):
     @pytest.mark.component
     def test_numerical_issues(self, system_frame):
         dt = DiagnosticsToolbox(system_frame)
-        warnings, next_steps = dt._collect_numerical_warnings()
-        assert len(warnings) == 3
+        warnings, _ = dt._collect_numerical_warnings()
+        assert len(warnings) == 1
         assert "WARNING: 3 Variables at or outside bounds (tol=0.0E+00)" in warnings
-        assert (
-            "WARNING: 2 Variables with extreme Jacobian values (<1.0E-08 or >1.0E+08)"
-            in warnings
-        )
-        assert (
-            "WARNING: 2 Constraints with extreme Jacobian values (<1.0E-08 or >1.0E+08)"
-            in warnings
-        )
 
     @pytest.mark.component
     def test_solve(self, system_frame):
@@ -201,17 +193,9 @@ def test_structural_issues(self, system_frame):
     @pytest.mark.component
     def test_numerical_issues(self, system_frame):
         dt = DiagnosticsToolbox(system_frame)
-        warnings, next_steps = dt._collect_numerical_warnings()
-        assert len(warnings) == 3
+        warnings, _ = dt._collect_numerical_warnings()
+        assert len(warnings) == 1
         assert "WARNING: 3 Variables at or outside bounds (tol=0.0E+00)" in warnings
-        assert (
-            "WARNING: 2 Variables with extreme Jacobian values (<1.0E-08 or >1.0E+08)"
-            in warnings
-        )
-        assert (
-            "WARNING: 2 Constraints with extreme Jacobian values (<1.0E-08 or >1.0E+08)"
-            in warnings
-        )
 
     @pytest.mark.component
     def test_solve(self, system_frame):
@@ -334,13 +318,30 @@ def test_numerical_issues(self, system_frame):
 
     @pytest.mark.solver
     @pytest.mark.component
-    def test_condition_number(self, system_frame):
+    @linux_platform_only
+    def test_condition_number_on_linux(self, system_frame):
         m = system_frame
 
         # Check condition number to confirm scaling
         jac, _ = get_jacobian(m, scaled=False)
-        assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            3.6954462e15, rel=1e-3
+        assert jacobian_cond(jac=jac, scaled=False) == pytest.approx(
+            2.987650e15, rel=1e-2
+        )
+
+    @pytest.mark.solver
+    @pytest.mark.component
+    @reference_platform_only
+    def test_condition_number_on_windows(self, system_frame):
+        m = system_frame
+
+        # Check condition number to confirm scaling
+        jac, _ = get_jacobian(m, scaled=False)
+        cond = jacobian_cond(jac=jac, scaled=False)
+        assert (
+            # Python 3.9-3.11
+            cond == pytest.approx(4.3021828e15, rel=1e-2)
+            # Python 3.12
+            or cond == pytest.approx(2.987651e15, rel=1e-2)
         )
 
 
@@ -377,7 +378,7 @@ def test_condition_number_on_linux(self, system_frame):
         jac, _ = get_jacobian(sm, scaled=False)
 
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            8.89313720973467e14, rel=1e-3
+            7.42017e14, rel=1e-3
         )
 
     @pytest.mark.solver
@@ -391,5 +392,5 @@ def test_condition_number_on_windows(self, system_frame):
         jac, _ = get_jacobian(sm, scaled=False)
 
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            8.82538591e14, rel=1e-2
+            7.43640e14, rel=1e-2
         )
diff --git a/watertap/flowsheets/full_water_resource_recovery_facility/test/test_full_WRRF_with_ASM1_ADM1.py b/watertap/flowsheets/full_water_resource_recovery_facility/test/test_full_WRRF_with_ASM1_ADM1.py
index 086f090518..49b307c1b0 100644
--- a/watertap/flowsheets/full_water_resource_recovery_facility/test/test_full_WRRF_with_ASM1_ADM1.py
+++ b/watertap/flowsheets/full_water_resource_recovery_facility/test/test_full_WRRF_with_ASM1_ADM1.py
@@ -26,6 +26,7 @@
 # Some more information about this module
 __author__ = "Alejandro Garciadiego, Xinhong Liu, Adam Atia, Marcus Holly"
 
+import platform
 import pytest
 
 from pyomo.environ import (
@@ -43,6 +44,27 @@
 
 import watertap.flowsheets.full_water_resource_recovery_facility.BSM2 as BSM2
 
+is_reference_platform = (
+    platform.system() == "Windows" and platform.python_version_tuple()[0] == "3"
+)
+is_linux_platform = (
+    platform.system() == "Linux" and platform.python_version_tuple()[0] == "3"
+)
+
+reference_platform_only = pytest.mark.xfail(
+    condition=(not is_reference_platform),
+    run=True,
+    strict=False,
+    reason="These tests are expected to pass only on the reference platform (Python 3 on Windows)",
+)
+
+linux_platform_only = pytest.mark.xfail(
+    condition=(not is_linux_platform),
+    run=True,
+    strict=False,
+    reason="These tests are expected to pass only on the Linux platform (Python 3)",
+)
+
 
 class TestFullFlowsheet:
     @pytest.fixture(scope="class")
@@ -182,7 +204,7 @@ def test_condition_number(self, system_frame):
         # Check condition number to confirm scaling
         jac, _ = get_jacobian(m.scaled_model, scaled=False)
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            5.4097699e9, rel=1e-3
+            6.80815e9, rel=1e-3
         )
 
     @pytest.mark.component
@@ -194,7 +216,28 @@ def test_display(self, system_frame):
 
     @pytest.mark.requires_idaes_solver
     @pytest.mark.component
-    def test_optimization(self, optimized_system_frame):
+    @reference_platform_only
+    def test_optimization_windows(self, optimized_system_frame):
+        m = optimized_system_frame
+        assert_optimal_termination(m.rescaled_results)
+
+        assert degrees_of_freedom(m) == 16
+
+        # Check condition number to confirm scaling
+        jac, _ = get_jacobian(m.rescaled_model, scaled=False)
+        cond = jacobian_cond(jac=jac, scaled=False)
+        assert (
+            # Python 3.9 and 3.10
+            cond == pytest.approx(1.95367e11, rel=1e-2)
+            # Python 3.11 and 3.12
+            # or cond == pytest.approx(3.44132e11, rel=1e-2)
+            or cond == pytest.approx(2.71713e11, rel=1e-2)
+        )
+
+    @pytest.mark.requires_idaes_solver
+    @pytest.mark.component
+    @linux_platform_only
+    def test_optimization_linux(self, optimized_system_frame):
         m = optimized_system_frame
         assert_optimal_termination(m.rescaled_results)
 
@@ -203,5 +246,7 @@ def test_optimization(self, optimized_system_frame):
         # Check condition number to confirm scaling
         jac, _ = get_jacobian(m.rescaled_model, scaled=False)
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            3.1695338e11, rel=1e-3
+            3.44152e11,
+            # 2.71713e11,
+            rel=1e-3,
         )
diff --git a/watertap/property_models/multicomp_aq_sol_prop_pack.py b/watertap/property_models/multicomp_aq_sol_prop_pack.py
index c8f76199e8..89e0465fe4 100644
--- a/watertap/property_models/multicomp_aq_sol_prop_pack.py
+++ b/watertap/property_models/multicomp_aq_sol_prop_pack.py
@@ -1407,8 +1407,8 @@ def _flow_mol_phase_comp(self):
 
             def rule_flow_mol_phase_comp(b, p, j):
                 return (
-                    b.flow_mass_phase_comp[p, j]
-                    == b.flow_mol_phase_comp[p, j] * b.params.mw_comp[j]
+                    b.flow_mass_phase_comp[p, j] / b.params.mw_comp[j]
+                    == b.flow_mol_phase_comp[p, j]
                 )
 
             self.eq_flow_mol_phase_comp = Constraint(
diff --git a/watertap/property_models/tests/test_multicomp_aq_sol_prop_pack.py b/watertap/property_models/tests/test_multicomp_aq_sol_prop_pack.py
index 8262b31a37..7e1b3c23d8 100644
--- a/watertap/property_models/tests/test_multicomp_aq_sol_prop_pack.py
+++ b/watertap/property_models/tests/test_multicomp_aq_sol_prop_pack.py
@@ -2667,6 +2667,7 @@ def test_scaler(self):
 
         assert len(m.fs.stream[0].scaling_factor) == 154
 
+        # Variables
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].temperature
         ] == pytest.approx(0.01, rel=1e-3)
@@ -2910,6 +2911,8 @@ def test_scaler(self):
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].total_hardness
         ] == pytest.approx(0.0001443, rel=1e-3)
+
+        # Constraints
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_conc_mol_phase_comp["Liq", "H2O"]
         ] == pytest.approx(0.001, rel=1e-3)
@@ -2954,19 +2957,19 @@ def test_scaler(self):
         ] == pytest.approx(1, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_mass_frac_phase_comp["Liq", "Ca_2+"]
-        ] == pytest.approx(2617, rel=1e-3)
+        ] == pytest.approx(1e4, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_mass_frac_phase_comp["Liq", "SO4_2-"]
-        ] == pytest.approx(468.1, rel=1e-3)
+        ] == pytest.approx(1e3, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_mass_frac_phase_comp["Liq", "Na_+"]
-        ] == pytest.approx(89.89, rel=1e-3)
+        ] == pytest.approx(1e2, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_mass_frac_phase_comp["Liq", "Cl_-"]
-        ] == pytest.approx(49.74, rel=1e-3)
+        ] == pytest.approx(1e2, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_mass_frac_phase_comp["Liq", "Mg_2+"]
-        ] == pytest.approx(717.2, rel=1e-3)
+        ] == pytest.approx(1e3, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_dens_mass_solvent
         ] == pytest.approx(0.001, rel=1e-3)
@@ -3005,28 +3008,28 @@ def test_scaler(self):
         ] == pytest.approx(18, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_flow_mol_phase_comp["Liq", "H2O"]
-        ] == pytest.approx(1.037, rel=1e-3)
+        ] == pytest.approx(0.018, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_flow_mol_phase_comp["Liq", "Ca_2+"]
-        ] == pytest.approx(2500, rel=1e-3)
+        ] == pytest.approx(100, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_flow_mol_phase_comp["Liq", "SO4_2-"]
-        ] == pytest.approx(468.2, rel=1e-3)
+        ] == pytest.approx(96, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_flow_mol_phase_comp["Liq", "Na_+"]
-        ] == pytest.approx(89.91, rel=1e-3)
+        ] == pytest.approx(2.3, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_flow_mol_phase_comp["Liq", "Cl_-"]
-        ] == pytest.approx(49.75, rel=1e-3)
+        ] == pytest.approx(3.5, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_flow_mol_phase_comp["Liq", "Mg_2+"]
-        ] == pytest.approx(717.4, rel=1e-3)
+        ] == pytest.approx(24, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_debye_huckel_constant
-        ] == pytest.approx(64.32, rel=1e-3)
+        ] == pytest.approx(12.49, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_enth_mass_phase["Liq"]
-        ] == pytest.approx(0.0000007982, rel=1e-3)
+        ] == pytest.approx(8.712e-07, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_total_dissolved_solids
         ] == pytest.approx(0.0000278, rel=1e-3)
@@ -3053,10 +3056,10 @@ def test_scaler(self):
         ] == pytest.approx(1000, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_pressure_osm_phase["Liq"]
-        ] == pytest.approx(0.0000009278, rel=1e-3)
+        ] == pytest.approx(1.290e-06, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_pressure_sat
-        ] == pytest.approx(0.0003216, rel=1e-3)
+        ] == pytest.approx(1e-3, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_visc_k_phase["Liq"]
         ] == pytest.approx(1000, rel=1e-3)
@@ -3095,19 +3098,19 @@ def test_scaler(self):
         ] == pytest.approx(1545, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_elec_mobility_phase_comp["Liq", "Ca_2+"]
-        ] == pytest.approx(16220000, rel=1e-3)
+        ] == pytest.approx(4.309e07, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_elec_mobility_phase_comp["Liq", "SO4_2-"]
-        ] == pytest.approx(12120000, rel=1e-3)
+        ] == pytest.approx(4.309e07, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_elec_mobility_phase_comp["Liq", "Na_+"]
-        ] == pytest.approx(19320000, rel=1e-3)
+        ] == pytest.approx(8.617e07, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_elec_mobility_phase_comp["Liq", "Cl_-"]
-        ] == pytest.approx(12660000, rel=1e-3)
+        ] == pytest.approx(8.617e07, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_elec_mobility_phase_comp["Liq", "Mg_2+"]
-        ] == pytest.approx(18200000, rel=1e-3)
+        ] == pytest.approx(4.309e07, rel=1e-3)
         assert m.fs.stream[0].scaling_factor[
             m.fs.stream[0].eq_elec_cond_phase["Liq"]
         ] == pytest.approx(3.225, rel=1e-3)
diff --git a/watertap/unit_models/reverse_osmosis_0D.py b/watertap/unit_models/reverse_osmosis_0D.py
index e5b8342075..479742d0a5 100644
--- a/watertap/unit_models/reverse_osmosis_0D.py
+++ b/watertap/unit_models/reverse_osmosis_0D.py
@@ -189,3 +189,15 @@ def calculate_scaling_factors(self):
         if hasattr(self, "width"):
             if iscale.get_scaling_factor(self.width) is None:
                 iscale.set_scaling_factor(self.width, 1)
+
+        for (t, x, j), condata in self.eq_mass_frac_permeate.items():
+            sf = iscale.get_scaling_factor(
+                self.flux_mass_phase_comp[t, x, "Liq", j], default=1
+            )
+            iscale.constraint_scaling_transform(condata, sf)
+
+        for (t, p, j), condata in self.eq_permeate_production.items():
+            sf = iscale.get_scaling_factor(
+                self.mixed_permeate[t].get_material_flow_terms(p, j), 1
+            )
+            iscale.constraint_scaling_transform(condata, sf)
diff --git a/watertap/unit_models/reverse_osmosis_1D.py b/watertap/unit_models/reverse_osmosis_1D.py
index 3840f173a3..1ab253c095 100644
--- a/watertap/unit_models/reverse_osmosis_1D.py
+++ b/watertap/unit_models/reverse_osmosis_1D.py
@@ -242,3 +242,59 @@ def calculate_scaling_factors(self):
             for v in self.deltaP.values():
                 if iscale.get_scaling_factor(v) is None:
                     iscale.set_scaling_factor(v, 1e-4)
+
+        if hasattr(self, "eq_permeate_outlet_isobaric"):
+            for (t, x), condata in self.eq_permeate_outlet_isobaric.items():
+                sf_P = min(
+                    iscale.get_scaling_factor(
+                        self.permeate_side[t, x].pressure, default=1
+                    ),
+                    iscale.get_scaling_factor(
+                        self.mixed_permeate[t].pressure, default=1
+                    ),
+                )
+                iscale.constraint_scaling_transform(condata, sf_P)
+
+        if self.config.pressure_change_type == PressureChangeType.fixed_per_stage:
+            # TODO additional scaling
+            pass
+        else:
+            for t, condata in self.eq_pressure_drop.items():
+                sf_dP = iscale.get_scaling_factor(self.deltaP[t])
+                iscale.constraint_scaling_transform(condata, sf_dP)
+
+        for (t, p, j), condata in self.eq_permeate_production.items():
+            sf = iscale.get_scaling_factor(
+                self.mixed_permeate[t].get_material_flow_terms(p, j), default=1
+            )
+            iscale.constraint_scaling_transform(condata, sf)
+
+        for (t, x, p, j), condata in self.eq_mass_transfer_term.items():
+            sf = min(
+                iscale.get_scaling_factor(
+                    self.mass_transfer_phase_comp[t, x, p, j], default=1
+                ),
+                iscale.get_scaling_factor(
+                    self.feed_side.mass_transfer_term[t, x, p, j], default=1
+                ),
+            )
+            iscale.constraint_scaling_transform(condata, sf)
+
+        for (t, x, p, j), condata in self.eq_connect_mass_transfer.items():
+            sf = iscale.get_scaling_factor(
+                self.permeate_side[t, x].get_material_flow_terms(p, j), default=1
+            )
+            iscale.constraint_scaling_transform(condata, sf)
+
+        for (t, x, p, j), condata in self.eq_mass_flux_equal_mass_transfer.items():
+            sf_xfer = iscale.get_scaling_factor(
+                self.feed_side.mass_transfer_term[t, x, p, j]
+            )
+            sf_ell = iscale.get_scaling_factor(
+                self.length, default=1 / value(self.length)
+            )
+            sf_A = iscale.get_scaling_factor(self.area, default=1 / value(self.area))
+            iscale.set_scaling_factor(
+                self.flux_mass_phase_comp[t, x, p, j], sf_xfer * sf_ell / sf_A
+            )
+            iscale.constraint_scaling_transform(condata, sf_xfer * sf_ell)
diff --git a/watertap/unit_models/reverse_osmosis_base.py b/watertap/unit_models/reverse_osmosis_base.py
index 91c176e68b..c34a9cdd4d 100644
--- a/watertap/unit_models/reverse_osmosis_base.py
+++ b/watertap/unit_models/reverse_osmosis_base.py
@@ -876,6 +876,9 @@ def calculate_scaling_factors(self):
                         )
                     )
                     iscale.set_scaling_factor(v, sf)
+                    iscale.constraint_scaling_transform(
+                        self.eq_flux_mass[t, x, p, j], sf
+                    )
                 elif comp.is_solute():  # scaling based on solute flux equation
                     sf = iscale.get_scaling_factor(
                         self.B_comp[t, j]
@@ -883,6 +886,15 @@ def calculate_scaling_factors(self):
                         self.feed_side.properties[t, x].conc_mass_phase_comp[p, j]
                     )
                     iscale.set_scaling_factor(v, sf)
+                    iscale.constraint_scaling_transform(
+                        self.eq_flux_mass[t, x, p, j], sf
+                    )
+
+        for (t, j), condata in self.eq_recovery_mass_phase_comp.items():
+            sf = iscale.get_scaling_factor(
+                self.mixed_permeate[t].flow_mass_phase_comp["Liq", j], default=1
+            )
+            iscale.constraint_scaling_transform(condata, sf)
 
     @property
     def default_costing_method(self):
diff --git a/watertap/unit_models/tests/test_aeration_tank.py b/watertap/unit_models/tests/test_aeration_tank.py
index be35626988..46b1ba5b4d 100644
--- a/watertap/unit_models/tests/test_aeration_tank.py
+++ b/watertap/unit_models/tests/test_aeration_tank.py
@@ -14,7 +14,6 @@
 Authors: Andrew Lee, Vibhav Dabadghao
 """
 
-from io import StringIO
 import pytest
 from pyomo.environ import (
     ConcreteModel,
@@ -41,7 +40,6 @@
     get_jacobian,
     jacobian_cond,
 )
-from idaes.core.scaling.scaler_profiling import ScalingProfiler
 import idaes.core.util.scaling as iscale
 from idaes.core.scaling.scaling_base import ScalerBase
 from idaes.core.util.testing import (
@@ -429,9 +427,9 @@ def test_constraint_scaling_routine(self, model):
 
         scaler.constraint_scaling_routine(model.fs.unit)
 
-        sfx_out = model.fs.unit.control_volume.properties_out[0].scaling_factor
-        assert isinstance(sfx_out, Suffix)
-        assert len(sfx_out) == 0
+        assert not hasattr(
+            model.fs.unit.control_volume.properties_out[0], "scaling_factor"
+        )
 
         sfx_rxn = model.fs.unit.control_volume.reactions[0].scaling_factor
         assert isinstance(sfx_rxn, Suffix)
@@ -641,7 +639,7 @@ def test_example_case_scaler(self):
         sm = TransformationFactory("core.scale_model").create_using(m, rename=False)
         jac, _ = get_jacobian(sm, scaled=False)
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            6.732817e6, rel=1e-3
+            1.1526931e7, rel=1e-3
         )
 
 
@@ -771,71 +769,72 @@ def perturb_solution(m):
     m.fs.unit.injection[0, "Liq", "S_O"].fix(2e-3 * 0.7)
 
 
-@pytest.mark.requires_idaes_solver
-@pytest.mark.unit
-def test_scaling_profiler_with_scalers():
-    sp = ScalingProfiler(
-        build_model=build_model,
-        user_scaling=scale_vars_with_scalers,
-        perturb_state=perturb_solution,
-    )
-
-    stream = StringIO()
-
-    sp.report_scaling_profiles(stream=stream)
-
-    expected = """
-============================================================================
-Scaling Profile Report
-----------------------------------------------------------------------------
-Scaling Method           || User Scaling           || Perfect Scaling
-Unscaled                 || 1.826E+16 | Solved 14  ||
-Vars Only                || 2.740E+17 | Solved 4   || 2.014E+21 | Solved 4  
-Harmonic                 || 2.740E+17 | Solved 4   || 4.443E+22 | Solved 27 
-Inverse Sum              || 2.740E+17 | Solved 4   || 2.399E+14 | Solved 4  
-Inverse Root Sum Squares || 2.740E+17 | Solved 4   || 3.412E+14 | Solved 4  
-Inverse Maximum          || 2.740E+17 | Solved 4   || 4.809E+14 | Solved 4  
-Inverse Minimum          || 2.740E+17 | Solved 4   || 4.455E+22 | Solved 24 
-Nominal L1 Norm          || 2.740E+17 | Solved 4   || 2.841E+14 | Solved 3  
-Nominal L2 Norm          || 2.740E+17 | Solved 4   || 3.755E+14 | Solved 3  
-Actual L1 Norm           || 2.740E+17 | Solved 4   || 5.461E+13 | Solved 4  
-Actual L2 Norm           || 2.740E+17 | Solved 4   || 6.491E+13 | Solved 4  
-============================================================================
-"""
-
-    assert stream.getvalue() == expected
-
-
-@pytest.mark.unit
-@pytest.mark.requires_idaes_solver
-def test_scaling_profiler_with_iscale():
-    sp = ScalingProfiler(
-        build_model=build_model,
-        user_scaling=scale_vars_with_iscale,
-        perturb_state=perturb_solution,
-    )
-
-    stream = StringIO()
-
-    sp.report_scaling_profiles(stream=stream)
-
-    expected = """
-============================================================================
-Scaling Profile Report
-----------------------------------------------------------------------------
-Scaling Method           || User Scaling           || Perfect Scaling
-Unscaled                 || 1.826E+16 | Solved 14  ||
-Vars Only                || 8.948E+12 | Solved 4   || 2.014E+21 | Solved 4  
-Harmonic                 || 1.044E+17 | Solved 44  || 4.443E+22 | Solved 27 
-Inverse Sum              || 5.247E+17 | Solved 66  || 2.399E+14 | Solved 4  
-Inverse Root Sum Squares || 5.220E+17 | Solved 73  || 3.412E+14 | Solved 4  
-Inverse Maximum          || 5.208E+17 | Solved 66  || 4.809E+14 | Solved 4  
-Inverse Minimum          || 2.103E+17 | Solved 85  || 4.455E+22 | Solved 24 
-Nominal L1 Norm          || 7.817E+09 | Solved 6   || 2.841E+14 | Solved 3  
-Nominal L2 Norm          || 1.278E+10 | Solved 6   || 3.755E+14 | Solved 3  
-Actual L1 Norm           || 3.950E+09 | Solved 3   || 5.461E+13 | Solved 4  
-Actual L2 Norm           || 4.339E+09 | Solved 3   || 6.491E+13 | Solved 4  
-============================================================================
-"""
-
-    assert stream.getvalue() == expected
+# TODO Replace these scaling profiler tests with more detailed convergence analysis
+# @pytest.mark.requires_idaes_solver
+# @pytest.mark.unit
+# def test_scaling_profiler_with_scalers():
+#     sp = ScalingProfiler(
+#         build_model=build_model,
+#         user_scaling=scale_vars_with_scalers,
+#         perturb_state=perturb_solution,
+#     )
+
+#     stream = StringIO()
+
+#     sp.report_scaling_profiles(stream=stream)
+
+#     expected = """
+# ============================================================================
+# Scaling Profile Report
+# ----------------------------------------------------------------------------
+# Scaling Method           || User Scaling           || Perfect Scaling
+# Unscaled                 || 1.826E+16 | Solved 14  ||
+# Vars Only                || 2.740E+17 | Solved 4   || 2.014E+21 | Solved 4
+# Harmonic                 || 2.740E+17 | Solved 4   || 4.443E+22 | Solved 27
+# Inverse Sum              || 2.740E+17 | Solved 4   || 2.399E+14 | Solved 4
+# Inverse Root Sum Squares || 2.740E+17 | Solved 4   || 3.412E+14 | Solved 4
+# Inverse Maximum          || 2.740E+17 | Solved 4   || 4.809E+14 | Solved 4
+# Inverse Minimum          || 2.740E+17 | Solved 4   || 4.455E+22 | Solved 24
+# Nominal L1 Norm          || 2.740E+17 | Solved 4   || 2.841E+14 | Solved 3
+# Nominal L2 Norm          || 2.740E+17 | Solved 4   || 3.755E+14 | Solved 3
+# Actual L1 Norm           || 2.740E+17 | Solved 4   || 5.461E+13 | Solved 4
+# Actual L2 Norm           || 2.740E+17 | Solved 4   || 6.491E+13 | Solved 4
+# ============================================================================
+# """
+
+#     assert stream.getvalue() == expected
+
+
+# @pytest.mark.unit
+# @pytest.mark.requires_idaes_solver
+# def test_scaling_profiler_with_iscale():
+#     sp = ScalingProfiler(
+#         build_model=build_model,
+#         user_scaling=scale_vars_with_iscale,
+#         perturb_state=perturb_solution,
+#     )
+
+#     stream = StringIO()
+
+#     sp.report_scaling_profiles(stream=stream)
+
+#     expected = """
+# ============================================================================
+# Scaling Profile Report
+# ----------------------------------------------------------------------------
+# Scaling Method           || User Scaling           || Perfect Scaling
+# Unscaled                 || 1.826E+16 | Solved 14  ||
+# Vars Only                || 8.948E+12 | Solved 4   || 2.014E+21 | Solved 4
+# Harmonic                 || 1.044E+17 | Solved 44  || 4.443E+22 | Solved 27
+# Inverse Sum              || 5.247E+17 | Solved 66  || 2.399E+14 | Solved 4
+# Inverse Root Sum Squares || 5.220E+17 | Solved 73  || 3.412E+14 | Solved 4
+# Inverse Maximum          || 5.208E+17 | Solved 66  || 4.809E+14 | Solved 4
+# Inverse Minimum          || 2.103E+17 | Solved 85  || 4.455E+22 | Solved 24
+# Nominal L1 Norm          || 7.817E+09 | Solved 6   || 2.841E+14 | Solved 3
+# Nominal L2 Norm          || 1.278E+10 | Solved 6   || 3.755E+14 | Solved 3
+# Actual L1 Norm           || 3.950E+09 | Solved 3   || 5.461E+13 | Solved 4
+# Actual L2 Norm           || 4.339E+09 | Solved 3   || 6.491E+13 | Solved 4
+# ============================================================================
+# """
+
+#     assert stream.getvalue() == expected
diff --git a/watertap/unit_models/tests/test_anaerobic_digester.py b/watertap/unit_models/tests/test_anaerobic_digester.py
index cb9b28bc92..54ea6d2176 100644
--- a/watertap/unit_models/tests/test_anaerobic_digester.py
+++ b/watertap/unit_models/tests/test_anaerobic_digester.py
@@ -351,9 +351,9 @@ def test_constraint_scaling_routine(self, model):
 
         scaler.constraint_scaling_routine(model.fs.unit)
 
-        sfx_out = model.fs.unit.liquid_phase.properties_out[0].scaling_factor
-        assert isinstance(sfx_out, Suffix)
-        assert len(sfx_out) == 0
+        assert not hasattr(
+            model.fs.unit.liquid_phase.properties_out[0], "scaling_factor"
+        )
 
         sfx_rxn = model.fs.unit.liquid_phase.reactions[0].scaling_factor
         assert isinstance(sfx_rxn, Suffix)
@@ -525,7 +525,7 @@ def test_example_case_scaler_scaling_default(self):
         sm = TransformationFactory("core.scale_model").create_using(m, rename=False)
         jac, _ = get_jacobian(sm, scaled=False)
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            3.41888877e11, rel=1e-3
+            2.504226e11, rel=1e-3
         )
 
     @pytest.mark.integration
@@ -610,5 +610,5 @@ def test_example_case_scaler_scaling(self):
         sm = TransformationFactory("core.scale_model").create_using(m, rename=False)
         jac, _ = get_jacobian(sm, scaled=False)
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            7.986901e10, rel=1e-3
+            8.432989e10, rel=1e-3
         )
diff --git a/watertap/unit_models/tests/test_clarifier.py b/watertap/unit_models/tests/test_clarifier.py
index 0f4a8aaf7d..e4c9f821a4 100644
--- a/watertap/unit_models/tests/test_clarifier.py
+++ b/watertap/unit_models/tests/test_clarifier.py
@@ -13,7 +13,6 @@
 Tests for clarifier.
 """
 __author__ = "Chenyu Wang"
-from io import StringIO
 import pytest
 from pyomo.environ import (
     ConcreteModel,
@@ -29,7 +28,6 @@
     get_jacobian,
     jacobian_cond,
 )
-from idaes.core.scaling.scaler_profiling import ScalingProfiler
 from watertap.core.solvers import get_solver
 
 from watertap.unit_models.tests.unit_test_harness import UnitTestHarness
@@ -570,7 +568,7 @@ def test_example_case_scaler(self):
         sm = TransformationFactory("core.scale_model").create_using(m, rename=False)
         jac, _ = get_jacobian(sm, scaled=False)
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            2.0028333e4, rel=1e-3
+            2729.85, rel=1e-3
         )
 
 
@@ -666,71 +664,72 @@ def perturb_solution(m):
     m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(27 * 1.5 * units.g / units.m**3)
 
 
-@pytest.mark.requires_idaes_solver
-@pytest.mark.unit
-def test_scaling_profiler_with_scalers():
-    sp = ScalingProfiler(
-        build_model=build_model,
-        user_scaling=scale_vars_with_scalers,
-        perturb_state=perturb_solution,
-    )
-
-    stream = StringIO()
-
-    sp.report_scaling_profiles(stream=stream)
-
-    expected = """
-============================================================================
-Scaling Profile Report
-----------------------------------------------------------------------------
-Scaling Method           || User Scaling           || Perfect Scaling
-Unscaled                 || 6.241E+06 | Solved 3   ||
-Vars Only                || 1.359E+10 | Solved 3   || 1.356E+14 | Solved 1  
-Harmonic                 || 1.359E+10 | Solved 3   || 1.228E+02 | Solved 1  
-Inverse Sum              || 1.359E+10 | Solved 3   || 6.636E+03 | Solved 1  
-Inverse Root Sum Squares || 1.359E+10 | Solved 3   || 6.633E+03 | Solved 1  
-Inverse Maximum          || 1.359E+10 | Solved 3   || 6.636E+03 | Solved 1  
-Inverse Minimum          || 1.359E+10 | Solved 3   || 9.327E+01 | Solved 1  
-Nominal L1 Norm          || 1.359E+10 | Solved 3   || 1.817E+03 | Solved 1  
-Nominal L2 Norm          || 1.359E+10 | Solved 3   || 3.354E+03 | Solved 1  
-Actual L1 Norm           || 1.359E+10 | Solved 3   || 9.450E+01 | Solved 1  
-Actual L2 Norm           || 1.359E+10 | Solved 3   || 8.480E+01 | Solved 1  
-============================================================================
-"""
-
-    assert stream.getvalue() == expected
-
-
-@pytest.mark.requires_idaes_solver
-@pytest.mark.unit
-def test_scaling_profiler_with_iscale():
-    sp = ScalingProfiler(
-        build_model=build_model,
-        user_scaling=scale_vars_with_iscale,
-        perturb_state=perturb_solution,
-    )
-
-    stream = StringIO()
-
-    sp.report_scaling_profiles(stream=stream)
-
-    expected = """
-============================================================================
-Scaling Profile Report
-----------------------------------------------------------------------------
-Scaling Method           || User Scaling           || Perfect Scaling
-Unscaled                 || 6.241E+06 | Solved 3   ||
-Vars Only                || 2.999E+09 | Solved 2   || 1.356E+14 | Solved 1  
-Harmonic                 || 1.027E+06 | Solved 2   || 1.228E+02 | Solved 1  
-Inverse Sum              || 8.148E+06 | Solved 2   || 6.636E+03 | Solved 1  
-Inverse Root Sum Squares || 8.114E+06 | Solved 2   || 6.633E+03 | Solved 1  
-Inverse Maximum          || 8.139E+06 | Solved 2   || 6.636E+03 | Solved 1  
-Inverse Minimum          || 1.257E+06 | Solved 2   || 9.327E+01 | Solved 1  
-Nominal L1 Norm          || 1.632E+09 | Solved 2   || 1.817E+03 | Solved 1  
-Nominal L2 Norm          || 1.972E+09 | Solved 2   || 3.354E+03 | Solved 1  
-Actual L1 Norm           || 1.776E+05 | Solved 2   || 9.450E+01 | Solved 1  
-Actual L2 Norm           || 1.723E+05 | Solved 2   || 8.480E+01 | Solved 1  
-============================================================================
-"""
-
-    assert stream.getvalue() == expected
+# TODO Replace these scaling profiler tests with more detailed convergence analysis
+# @pytest.mark.requires_idaes_solver
+# @pytest.mark.unit
+# def test_scaling_profiler_with_scalers():
+#     sp = ScalingProfiler(
+#         build_model=build_model,
+#         user_scaling=scale_vars_with_scalers,
+#         perturb_state=perturb_solution,
+#     )
+
+#     stream = StringIO()
+
+#     sp.report_scaling_profiles(stream=stream)
+
+#     expected = """
+# ============================================================================
+# Scaling Profile Report
+# ----------------------------------------------------------------------------
+# Scaling Method           || User Scaling           || Perfect Scaling
+# Unscaled                 || 6.241E+06 | Solved 3   ||
+# Vars Only                || 1.359E+10 | Solved 3   || 1.356E+14 | Solved 1
+# Harmonic                 || 1.359E+10 | Solved 3   || 1.228E+02 | Solved 1
+# Inverse Sum              || 1.359E+10 | Solved 3   || 6.636E+03 | Solved 1
+# Inverse Root Sum Squares || 1.359E+10 | Solved 3   || 6.633E+03 | Solved 1
+# Inverse Maximum          || 1.359E+10 | Solved 3   || 6.636E+03 | Solved 1
+# Inverse Minimum          || 1.359E+10 | Solved 3   || 9.327E+01 | Solved 1
+# Nominal L1 Norm          || 1.359E+10 | Solved 3   || 1.817E+03 | Solved 1
+# Nominal L2 Norm          || 1.359E+10 | Solved 3   || 3.354E+03 | Solved 1
+# Actual L1 Norm           || 1.359E+10 | Solved 3   || 9.450E+01 | Solved 1
+# Actual L2 Norm           || 1.359E+10 | Solved 3   || 8.480E+01 | Solved 1
+# ============================================================================
+# """
+
+#     assert stream.getvalue() == expected
+
+
+# @pytest.mark.requires_idaes_solver
+# @pytest.mark.unit
+# def test_scaling_profiler_with_iscale():
+#     sp = ScalingProfiler(
+#         build_model=build_model,
+#         user_scaling=scale_vars_with_iscale,
+#         perturb_state=perturb_solution,
+#     )
+
+#     stream = StringIO()
+
+#     sp.report_scaling_profiles(stream=stream)
+
+#     expected = """
+# ============================================================================
+# Scaling Profile Report
+# ----------------------------------------------------------------------------
+# Scaling Method           || User Scaling           || Perfect Scaling
+# Unscaled                 || 6.241E+06 | Solved 3   ||
+# Vars Only                || 2.999E+09 | Solved 2   || 1.356E+14 | Solved 1
+# Harmonic                 || 1.027E+06 | Solved 2   || 1.228E+02 | Solved 1
+# Inverse Sum              || 8.148E+06 | Solved 2   || 6.636E+03 | Solved 1
+# Inverse Root Sum Squares || 8.114E+06 | Solved 2   || 6.633E+03 | Solved 1
+# Inverse Maximum          || 8.139E+06 | Solved 2   || 6.636E+03 | Solved 1
+# Inverse Minimum          || 1.257E+06 | Solved 2   || 9.327E+01 | Solved 1
+# Nominal L1 Norm          || 1.632E+09 | Solved 2   || 1.817E+03 | Solved 1
+# Nominal L2 Norm          || 1.972E+09 | Solved 2   || 3.354E+03 | Solved 1
+# Actual L1 Norm           || 1.776E+05 | Solved 2   || 9.450E+01 | Solved 1
+# Actual L2 Norm           || 1.723E+05 | Solved 2   || 8.480E+01 | Solved 1
+# ============================================================================
+# """
+
+#     assert stream.getvalue() == expected
diff --git a/watertap/unit_models/tests/test_cstr.py b/watertap/unit_models/tests/test_cstr.py
index 9888b393db..63f54427de 100644
--- a/watertap/unit_models/tests/test_cstr.py
+++ b/watertap/unit_models/tests/test_cstr.py
@@ -13,7 +13,6 @@
 Tests for CSTR unit model.
 Authors: Marcus Holly
 """
-from io import StringIO
 import pytest
 
 from pyomo.environ import (
@@ -35,7 +34,6 @@
     get_jacobian,
     jacobian_cond,
 )
-from idaes.core.scaling.scaler_profiling import ScalingProfiler
 from watertap.unit_models.cstr import CSTR, CSTRScaler
 from watertap.costing import WaterTAPCosting
 
@@ -486,9 +484,9 @@ def test_constraint_scaling_routine(self, model):
 
         scaler.constraint_scaling_routine(model.fs.unit)
 
-        sfx_out = model.fs.unit.control_volume.properties_out[0].scaling_factor
-        assert isinstance(sfx_out, Suffix)
-        assert len(sfx_out) == 0
+        assert not hasattr(
+            model.fs.unit.control_volume.properties_out[0], "scaling_factor"
+        )
 
         sfx_rxn = model.fs.unit.control_volume.reactions[0].scaling_factor
         assert isinstance(sfx_rxn, Suffix)
@@ -689,7 +687,7 @@ def test_example_case_scaler(self):
         sm = TransformationFactory("core.scale_model").create_using(m, rename=False)
         jac, _ = get_jacobian(sm, scaled=False)
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            1.28432e5, rel=1e-3
+            1.65862e5, rel=1e-3
         )
 
 
@@ -810,71 +808,72 @@ def perturb_solution(m):
     m.fs.unit.volume[0].fix(1000 * 0.75 * units.m**3)
 
 
-@pytest.mark.requires_idaes_solver
-@pytest.mark.unit
-def test_scaling_profiler_with_scalers():
-    sp = ScalingProfiler(
-        build_model=build_model,
-        user_scaling=scale_vars_with_scalers,
-        perturb_state=perturb_solution,
-    )
-
-    stream = StringIO()
-
-    sp.report_scaling_profiles(stream=stream)
-
-    expected = """
-============================================================================
-Scaling Profile Report
-----------------------------------------------------------------------------
-Scaling Method           || User Scaling           || Perfect Scaling
-Unscaled                 || 1.196E+12 | Solved 4   ||
-Vars Only                || 1.801E+07 | Solved 5   || 1.674E+17 | Solved 1  
-Harmonic                 || 1.801E+07 | Solved 5   || 1.161E+05 | Solved 3  
-Inverse Sum              || 1.801E+07 | Solved 5   || 2.054E+02 | Solved 3  
-Inverse Root Sum Squares || 1.801E+07 | Solved 5   || 2.131E+02 | Solved 3  
-Inverse Maximum          || 1.801E+07 | Solved 5   || 2.427E+02 | Solved 3  
-Inverse Minimum          || 1.801E+07 | Solved 5   || 1.317E+05 | Solved 3  
-Nominal L1 Norm          || 1.801E+07 | Solved 5   || 2.185E+02 | Solved 3  
-Nominal L2 Norm          || 1.801E+07 | Solved 5   || 2.053E+02 | Solved 3  
-Actual L1 Norm           || 1.801E+07 | Solved 5   || 2.146E+02 | Solved 3  
-Actual L2 Norm           || 1.801E+07 | Solved 5   || 2.114E+02 | Solved 3  
-============================================================================
-"""
-
-    assert stream.getvalue() == expected
-
-
-@pytest.mark.requires_idaes_solver
-@pytest.mark.unit
-def test_scaling_profiler_with_iscale():
-    sp = ScalingProfiler(
-        build_model=build_model,
-        user_scaling=scale_vars_with_iscale,
-        perturb_state=perturb_solution,
-    )
-
-    stream = StringIO()
-
-    sp.report_scaling_profiles(stream=stream)
-
-    expected = """
-============================================================================
-Scaling Profile Report
-----------------------------------------------------------------------------
-Scaling Method           || User Scaling           || Perfect Scaling
-Unscaled                 || 1.196E+12 | Solved 4   ||
-Vars Only                || 3.003E+13 | Solved 3   || 1.674E+17 | Solved 1  
-Harmonic                 || 9.790E+09 | Solved 5   || 1.161E+05 | Solved 3  
-Inverse Sum              || 2.885E+07 | Solved 5   || 2.054E+02 | Solved 3  
-Inverse Root Sum Squares || 3.046E+07 | Solved 5   || 2.131E+02 | Solved 3  
-Inverse Maximum          || 3.164E+07 | Solved 5   || 2.427E+02 | Solved 3  
-Inverse Minimum          || 1.879E+10 | Solved 5   || 1.317E+05 | Solved 3  
-Nominal L1 Norm          || 1.261E+09 | Solved 5   || 2.185E+02 | Solved 3  
-Nominal L2 Norm          || 9.543E+08 | Solved 5   || 2.053E+02 | Solved 3  
-Actual L1 Norm           || 1.665E+06 | Solved 4   || 2.146E+02 | Solved 3  
-Actual L2 Norm           || 1.736E+06 | Solved 4   || 2.114E+02 | Solved 3  
-============================================================================
-"""
-
-    assert stream.getvalue() == expected
+# TODO Replace these scaling profiler tests with more detailed convergence analysis
+# @pytest.mark.requires_idaes_solver
+# @pytest.mark.unit
+# def test_scaling_profiler_with_scalers():
+#     sp = ScalingProfiler(
+#         build_model=build_model,
+#         user_scaling=scale_vars_with_scalers,
+#         perturb_state=perturb_solution,
+#     )
+
+#     stream = StringIO()
+
+#     sp.report_scaling_profiles(stream=stream)
+
+#     expected = """
+# ============================================================================
+# Scaling Profile Report
+# ----------------------------------------------------------------------------
+# Scaling Method           || User Scaling           || Perfect Scaling
+# Unscaled                 || 1.196E+12 | Solved 4   ||
+# Vars Only                || 1.801E+07 | Solved 5   || 1.674E+17 | Solved 1
+# Harmonic                 || 1.801E+07 | Solved 5   || 1.161E+05 | Solved 3
+# Inverse Sum              || 1.801E+07 | Solved 5   || 2.054E+02 | Solved 3
+# Inverse Root Sum Squares || 1.801E+07 | Solved 5   || 2.131E+02 | Solved 3
+# Inverse Maximum          || 1.801E+07 | Solved 5   || 2.427E+02 | Solved 3
+# Inverse Minimum          || 1.801E+07 | Solved 5   || 1.317E+05 | Solved 3
+# Nominal L1 Norm          || 1.801E+07 | Solved 5   || 2.185E+02 | Solved 3
+# Nominal L2 Norm          || 1.801E+07 | Solved 5   || 2.053E+02 | Solved 3
+# Actual L1 Norm           || 1.801E+07 | Solved 5   || 2.146E+02 | Solved 3
+# Actual L2 Norm           || 1.801E+07 | Solved 5   || 2.114E+02 | Solved 3
+# ============================================================================
+# """
+
+#     assert stream.getvalue() == expected
+
+
+# @pytest.mark.requires_idaes_solver
+# @pytest.mark.unit
+# def test_scaling_profiler_with_iscale():
+#     sp = ScalingProfiler(
+#         build_model=build_model,
+#         user_scaling=scale_vars_with_iscale,
+#         perturb_state=perturb_solution,
+#     )
+
+#     stream = StringIO()
+
+#     sp.report_scaling_profiles(stream=stream)
+
+#     expected = """
+# ============================================================================
+# Scaling Profile Report
+# ----------------------------------------------------------------------------
+# Scaling Method           || User Scaling           || Perfect Scaling
+# Unscaled                 || 1.196E+12 | Solved 4   ||
+# Vars Only                || 3.003E+13 | Solved 3   || 1.674E+17 | Solved 1
+# Harmonic                 || 9.790E+09 | Solved 5   || 1.161E+05 | Solved 3
+# Inverse Sum              || 2.885E+07 | Solved 5   || 2.054E+02 | Solved 3
+# Inverse Root Sum Squares || 3.046E+07 | Solved 5   || 2.131E+02 | Solved 3
+# Inverse Maximum          || 3.164E+07 | Solved 5   || 2.427E+02 | Solved 3
+# Inverse Minimum          || 1.879E+10 | Solved 5   || 1.317E+05 | Solved 3
+# Nominal L1 Norm          || 1.261E+09 | Solved 5   || 2.185E+02 | Solved 3
+# Nominal L2 Norm          || 9.543E+08 | Solved 5   || 2.053E+02 | Solved 3
+# Actual L1 Norm           || 1.665E+06 | Solved 4   || 2.146E+02 | Solved 3
+# Actual L2 Norm           || 1.736E+06 | Solved 4   || 2.114E+02 | Solved 3
+# ============================================================================
+# """
+
+#     assert stream.getvalue() == expected
diff --git a/watertap/unit_models/tests/test_cstr_injection.py b/watertap/unit_models/tests/test_cstr_injection.py
index 14ea2c05f1..9f0bca2b46 100644
--- a/watertap/unit_models/tests/test_cstr_injection.py
+++ b/watertap/unit_models/tests/test_cstr_injection.py
@@ -14,7 +14,6 @@
 Authors: Andrew Lee, Adam Atia, Vibhav Dabadghao
 """
 
-from io import StringIO
 import pytest
 from pyomo.environ import (
     ConcreteModel,
@@ -34,7 +33,6 @@
     get_jacobian,
     jacobian_cond,
 )
-from idaes.core.scaling.scaler_profiling import ScalingProfiler
 from idaes.core.scaling.scaling_base import ScalerBase
 from idaes.models.properties.examples.saponification_thermo import (
     SaponificationParameterBlock,
@@ -563,9 +561,9 @@ def test_constraint_scaling_routine(self, model):
 
         scaler.constraint_scaling_routine(model.fs.unit)
 
-        sfx_out = model.fs.unit.control_volume.properties_out[0].scaling_factor
-        assert isinstance(sfx_out, Suffix)
-        assert len(sfx_out) == 0
+        assert not hasattr(
+            model.fs.unit.control_volume.properties_out[0], "scaling_factor"
+        )
 
         sfx_rxn = model.fs.unit.control_volume.reactions[0].scaling_factor
         assert isinstance(sfx_rxn, Suffix)
@@ -777,7 +775,7 @@ def test_example_case_scaler(self):
         sm = TransformationFactory("core.scale_model").create_using(m, rename=False)
         jac, _ = get_jacobian(sm, scaled=False)
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            6.732817e6, rel=1e-3
+            1.1526931e7, rel=1e-3
         )
 
 
@@ -907,71 +905,72 @@ def perturb_solution(m):
     m.fs.unit.volume.fix(500 * 0.85)
 
 
-@pytest.mark.requires_idaes_solver
-@pytest.mark.unit
-def test_scaling_profiler_with_scalers():
-    sp = ScalingProfiler(
-        build_model=build_model,
-        user_scaling=scale_vars_with_scalers,
-        perturb_state=perturb_solution,
-    )
-
-    stream = StringIO()
-
-    sp.report_scaling_profiles(stream=stream)
-
-    expected = """
-============================================================================
-Scaling Profile Report
-----------------------------------------------------------------------------
-Scaling Method           || User Scaling           || Perfect Scaling
-Unscaled                 || 1.826E+16 | Solved 4   ||
-Vars Only                || 2.740E+17 | Solved 4   || 2.014E+21 | Solved 4  
-Harmonic                 || 2.740E+17 | Solved 4   || 4.443E+22 | Solved 18 
-Inverse Sum              || 2.740E+17 | Solved 4   || 2.399E+14 | Solved 4  
-Inverse Root Sum Squares || 2.740E+17 | Solved 4   || 3.412E+14 | Solved 4  
-Inverse Maximum          || 2.740E+17 | Solved 4   || 4.809E+14 | Solved 4  
-Inverse Minimum          || 2.740E+17 | Solved 4   || 4.455E+22 | Solved 18 
-Nominal L1 Norm          || 2.740E+17 | Solved 4   || 2.841E+14 | Solved 4  
-Nominal L2 Norm          || 2.740E+17 | Solved 4   || 3.755E+14 | Solved 4  
-Actual L1 Norm           || 2.740E+17 | Solved 4   || 5.461E+13 | Solved 4  
-Actual L2 Norm           || 2.740E+17 | Solved 4   || 6.491E+13 | Solved 4  
-============================================================================
-"""
-
-    assert stream.getvalue() == expected
-
-
-@pytest.mark.requires_idaes_solver
-@pytest.mark.unit
-def test_scaling_profiler_with_iscale():
-    sp = ScalingProfiler(
-        build_model=build_model,
-        user_scaling=scale_vars_with_iscale,
-        perturb_state=perturb_solution,
-    )
-
-    stream = StringIO()
-
-    sp.report_scaling_profiles(stream=stream)
-
-    expected = """
-============================================================================
-Scaling Profile Report
-----------------------------------------------------------------------------
-Scaling Method           || User Scaling           || Perfect Scaling
-Unscaled                 || 1.826E+16 | Solved 4   ||
-Vars Only                || 8.948E+12 | Solved 4   || 2.014E+21 | Solved 4  
-Harmonic                 || 1.044E+17 | Solved 57  || 4.443E+22 | Solved 18 
-Inverse Sum              || 5.247E+17 | Failed 50  || 2.399E+14 | Solved 4  
-Inverse Root Sum Squares || 5.220E+17 | Failed 55  || 3.412E+14 | Solved 4  
-Inverse Maximum          || 5.208E+17 | Failed 52  || 4.809E+14 | Solved 4  
-Inverse Minimum          || 2.103E+17 | Solved 65  || 4.455E+22 | Solved 18 
-Nominal L1 Norm          || 7.817E+09 | Solved 4   || 2.841E+14 | Solved 4  
-Nominal L2 Norm          || 1.278E+10 | Solved 4   || 3.755E+14 | Solved 4  
-Actual L1 Norm           || 3.950E+09 | Solved 3   || 5.461E+13 | Solved 4  
-Actual L2 Norm           || 4.339E+09 | Solved 3   || 6.491E+13 | Solved 4  
-============================================================================
-"""
-
-    assert stream.getvalue() == expected
+# TODO Replace these scaling profiler tests with more detailed convergence analysis
+# @pytest.mark.requires_idaes_solver
+# @pytest.mark.unit
+# def test_scaling_profiler_with_scalers():
+#     sp = ScalingProfiler(
+#         build_model=build_model,
+#         user_scaling=scale_vars_with_scalers,
+#         perturb_state=perturb_solution,
+#     )
+
+#     stream = StringIO()
+
+#     sp.report_scaling_profiles(stream=stream)
+
+#     expected = """
+# ============================================================================
+# Scaling Profile Report
+# ----------------------------------------------------------------------------
+# Scaling Method           || User Scaling           || Perfect Scaling
+# Unscaled                 || 1.826E+16 | Solved 4   ||
+# Vars Only                || 2.740E+17 | Solved 4   || 2.014E+21 | Solved 4
+# Harmonic                 || 2.740E+17 | Solved 4   || 4.443E+22 | Solved 18
+# Inverse Sum              || 2.740E+17 | Solved 4   || 2.399E+14 | Solved 4
+# Inverse Root Sum Squares || 2.740E+17 | Solved 4   || 3.412E+14 | Solved 4
+# Inverse Maximum          || 2.740E+17 | Solved 4   || 4.809E+14 | Solved 4
+# Inverse Minimum          || 2.740E+17 | Solved 4   || 4.455E+22 | Solved 18
+# Nominal L1 Norm          || 2.740E+17 | Solved 4   || 2.841E+14 | Solved 4
+# Nominal L2 Norm          || 2.740E+17 | Solved 4   || 3.755E+14 | Solved 4
+# Actual L1 Norm           || 2.740E+17 | Solved 4   || 5.461E+13 | Solved 4
+# Actual L2 Norm           || 2.740E+17 | Solved 4   || 6.491E+13 | Solved 4
+# ============================================================================
+# """
+
+#     assert stream.getvalue() == expected
+
+
+# @pytest.mark.requires_idaes_solver
+# @pytest.mark.unit
+# def test_scaling_profiler_with_iscale():
+#     sp = ScalingProfiler(
+#         build_model=build_model,
+#         user_scaling=scale_vars_with_iscale,
+#         perturb_state=perturb_solution,
+#     )
+
+#     stream = StringIO()
+
+#     sp.report_scaling_profiles(stream=stream)
+
+#     expected = """
+# ============================================================================
+# Scaling Profile Report
+# ----------------------------------------------------------------------------
+# Scaling Method           || User Scaling           || Perfect Scaling
+# Unscaled                 || 1.826E+16 | Solved 4   ||
+# Vars Only                || 8.948E+12 | Solved 4   || 2.014E+21 | Solved 4
+# Harmonic                 || 1.044E+17 | Solved 57  || 4.443E+22 | Solved 18
+# Inverse Sum              || 5.247E+17 | Failed 50  || 2.399E+14 | Solved 4
+# Inverse Root Sum Squares || 5.220E+17 | Failed 55  || 3.412E+14 | Solved 4
+# Inverse Maximum          || 5.208E+17 | Failed 52  || 4.809E+14 | Solved 4
+# Inverse Minimum          || 2.103E+17 | Solved 65  || 4.455E+22 | Solved 18
+# Nominal L1 Norm          || 7.817E+09 | Solved 4   || 2.841E+14 | Solved 4
+# Nominal L2 Norm          || 1.278E+10 | Solved 4   || 3.755E+14 | Solved 4
+# Actual L1 Norm           || 3.950E+09 | Solved 3   || 5.461E+13 | Solved 4
+# Actual L2 Norm           || 4.339E+09 | Solved 3   || 6.491E+13 | Solved 4
+# ============================================================================
+# """
+
+#     assert stream.getvalue() == expected
diff --git a/watertap/unit_models/tests/test_dewatering_unit.py b/watertap/unit_models/tests/test_dewatering_unit.py
index 10014a23b2..f6c610a6fe 100644
--- a/watertap/unit_models/tests/test_dewatering_unit.py
+++ b/watertap/unit_models/tests/test_dewatering_unit.py
@@ -33,7 +33,6 @@
     get_jacobian,
     jacobian_cond,
 )
-from idaes.core.scaling.scaling_base import ScalerBase
 from idaes.core.scaling.scaler_profiling import ScalingProfiler
 
 from idaes.models.unit_models.separator import SplittingType
@@ -1023,8 +1022,6 @@ def test_example_case_scaler(self):
 
         m.fs.unit.hydraulic_retention_time.fix()
 
-        sb = ScalerBase()
-
         scaler = DewatererScaler()
         scaler.scale_model(
             m.fs.unit,
@@ -1039,7 +1036,7 @@ def test_example_case_scaler(self):
         sm = TransformationFactory("core.scale_model").create_using(m, rename=False)
         jac, _ = get_jacobian(sm, scaled=False)
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            2.10895296e4, rel=1e-3
+            53957.5009, rel=1e-3
         )
 
 
diff --git a/watertap/unit_models/tests/test_nanofiltration_0D.py b/watertap/unit_models/tests/test_nanofiltration_0D.py
index d7fa0b585b..daccebf33a 100644
--- a/watertap/unit_models/tests/test_nanofiltration_0D.py
+++ b/watertap/unit_models/tests/test_nanofiltration_0D.py
@@ -1102,26 +1102,26 @@ def test_scaling(self, mcas_case):
 
         assert mcas_case.fs.unit.scaling_factor[
             mcas_case.fs.unit.material_balances[0.0, "H2O"]
-        ] == pytest.approx(0.01866, rel=1e-3)
+        ] == pytest.approx(0.1, rel=1e-3)
         assert mcas_case.fs.unit.scaling_factor[
             mcas_case.fs.unit.material_balances[0.0, "Ca_2+"]
-        ] == pytest.approx(104.7, rel=1e-3)
+        ] == pytest.approx(1e3, rel=1e-3)
         assert mcas_case.fs.unit.scaling_factor[
             mcas_case.fs.unit.material_balances[0.0, "SO4_2-"]
-        ] == pytest.approx(44.94, rel=1e-3)
+        ] == pytest.approx(1e2, rel=1e-3)
         assert mcas_case.fs.unit.scaling_factor[
             mcas_case.fs.unit.material_balances[0.0, "Mg_2+"]
-        ] == pytest.approx(1.722, rel=1e-3)
+        ] == pytest.approx(10, rel=1e-3)
         assert mcas_case.fs.unit.scaling_factor[
             mcas_case.fs.unit.material_balances[0.0, "Na_+"]
-        ] == pytest.approx(2.068, rel=1e-3)
+        ] == pytest.approx(10, rel=1e-3)
         assert mcas_case.fs.unit.scaling_factor[
             mcas_case.fs.unit.material_balances[0.0, "Cl_-"]
-        ] == pytest.approx(0.6174, rel=1e-3)
+        ] == pytest.approx(1, rel=1e-3)
 
         assert mcas_case.fs.unit.scaling_factor[
             mcas_case.fs.unit.separation_constraint[0.0, "Liq", "H2O"]
-        ] == pytest.approx(0.0233193, rel=1e-5)
+        ] == pytest.approx(0.1, rel=1e-5)
         assert mcas_case.fs.unit.scaling_factor[
             mcas_case.fs.unit.separation_constraint[0.0, "Liq", "Ca_2+"]
         ] == pytest.approx(0.004, rel=1e-5)
@@ -1143,13 +1143,13 @@ def test_scaling(self, mcas_case):
         )
         assert mcas_case.fs.unit.scaling_factor[
             mcas_case.fs.unit.retentate_pressure_balance[0.0]
-        ] == pytest.approx(0.0000025, rel=1e-3)
+        ] == pytest.approx(1e-5, rel=1e-3)
         assert mcas_case.fs.unit.scaling_factor[
             mcas_case.fs.unit.retentate_temperature_equality[0.0]
-        ] == pytest.approx(0.003354, rel=1e-3)
+        ] == pytest.approx(1e-2, rel=1e-3)
         assert mcas_case.fs.unit.scaling_factor[
             mcas_case.fs.unit.permeate_temperature_equality[0.0]
-        ] == pytest.approx(0.003354, rel=1e-3)
+        ] == pytest.approx(1e-2, rel=1e-3)
 
     @pytest.mark.requires_idaes_solver
     @pytest.mark.component
diff --git a/watertap/unit_models/tests/test_reverse_osmosis_1D.py b/watertap/unit_models/tests/test_reverse_osmosis_1D.py
index 7f6ac62307..63b788e61b 100644
--- a/watertap/unit_models/tests/test_reverse_osmosis_1D.py
+++ b/watertap/unit_models/tests/test_reverse_osmosis_1D.py
@@ -154,6 +154,9 @@ def configure(self):
             },
         }
 
+        self.skip_badly_scaled_vars = True
+        self.condition_number = 3.9140627e8  # Was previously 4.8157167e15
+
         return m
 
 
@@ -263,6 +266,9 @@ def configure(self):
             },
         }
 
+        self.skip_badly_scaled_vars = True
+        self.condition_number = 4.434273e9  # Was previously 3.40497838e14
+
         return m
 
 
@@ -388,6 +394,9 @@ def configure(self):
             },
         }
 
+        self.skip_badly_scaled_vars = True
+        self.condition_number = 1.2297005e11  # Was previously 4.14262023e17
+
         return m
 
 
@@ -498,6 +507,9 @@ def configure(self):
             },
         }
 
+        self.skip_badly_scaled_vars = True
+        self.condition_number = 1.161779e10  # Was previously 1.9452857e15
+
         return m
 
 
@@ -608,6 +620,9 @@ def configure(self):
             },
         }
 
+        self.skip_badly_scaled_vars = True
+        self.condition_number = 5.781112e9  # Was previously 6.166070e14
+
         return m
 
 
@@ -718,6 +733,9 @@ def configure(self):
             },
         }
 
+        self.skip_badly_scaled_vars = True
+        self.condition_number = 7.540701e9  # Was previously 1.27214033e17
+
         return m
 
 
@@ -830,6 +848,9 @@ def configure(self):
             },
         }
 
+        self.skip_badly_scaled_vars = True
+        self.condition_number = 1.693107e8  # Was previously 2.7895630e14
+
         return m
 
 
@@ -944,6 +965,9 @@ def configure(self):
             },
         }
 
+        self.skip_badly_scaled_vars = True
+        self.condition_number = 8.2477371e9  # Was previously 1.3695196e17
+
         return m
 
 
@@ -1059,4 +1083,7 @@ def configure(self):
             },
         }
 
+        self.skip_badly_scaled_vars = True
+        self.condition_number = 7.502876 * 1e10  # Was previously 1.5101144e17
+
         return m
diff --git a/watertap/unit_models/tests/test_thickener_unit.py b/watertap/unit_models/tests/test_thickener_unit.py
index 9485cd171a..ff7eadedec 100644
--- a/watertap/unit_models/tests/test_thickener_unit.py
+++ b/watertap/unit_models/tests/test_thickener_unit.py
@@ -690,7 +690,7 @@ def test_example_case_scaler(self):
         sm = TransformationFactory("core.scale_model").create_using(m, rename=False)
         jac, _ = get_jacobian(sm, scaled=False)
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            4.83364614e2, rel=1e-3
+            4723.49, rel=1e-3
         )
 
 
diff --git a/watertap/unit_models/tests/unit_test_harness.py b/watertap/unit_models/tests/unit_test_harness.py
index 64eead42cc..0751bfa765 100644
--- a/watertap/unit_models/tests/unit_test_harness.py
+++ b/watertap/unit_models/tests/unit_test_harness.py
@@ -21,6 +21,7 @@
 from watertap.core.solvers import get_solver
 from idaes.core.util.testing import initialization_tester
 import idaes.core.util.scaling as iscale
+from idaes.core.scaling.util import jacobian_cond
 import idaes.logger as idaeslog
 
 
@@ -179,6 +180,23 @@ def test_unit_solutions(self, frame):
             )
         results = opt.solve(blk, tee=True)
 
+        if hasattr(self, "condition_number"):
+            cond = jacobian_cond(blk)
+            if not cond == pytest.approx(
+                self.condition_number,
+                abs=self.default_absolute_tolerance,
+                rel=self.default_relative_tolerance,
+            ):
+                # Testing to a relative tolerance of 1e-3, so four significant digits should suffice.
+                raise AssertionError(
+                    f"Condition number: Expected {self.condition_number:.4E}, got {cond:.4E} instead"
+                )
+
+        if hasattr(self, "skip_badly_scaled_vars"):
+            skip_badly_scaled_vars = self.skip_badly_scaled_vars
+        else:
+            skip_badly_scaled_vars = False
+
         # check solve
         badly_scaled_vars = list(
             iscale.badly_scaled_var_generator(
@@ -188,7 +206,7 @@ def test_unit_solutions(self, frame):
                 zero=self.default_zero,
             )
         )
-        if badly_scaled_vars:
+        if not skip_badly_scaled_vars and len(badly_scaled_vars) > 0:
             lines = [
                 f"{x[0].name}\t{x[0].value}\tsf: {iscale.get_scaling_factor(x[0])}"
                 for x in badly_scaled_vars
diff --git a/watertap/unit_models/translators/tests/test_translator_adm1_asm1.py b/watertap/unit_models/translators/tests/test_translator_adm1_asm1.py
index 61837ae519..092c35fd1d 100644
--- a/watertap/unit_models/translators/tests/test_translator_adm1_asm1.py
+++ b/watertap/unit_models/translators/tests/test_translator_adm1_asm1.py
@@ -566,5 +566,5 @@ def test_example_case_scaler(self):
         sm = TransformationFactory("core.scale_model").create_using(m, rename=False)
         jac, _ = get_jacobian(sm, scaled=False)
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            5.4739117e3, rel=1e-3
+            3468.97, rel=1e-3
         )
diff --git a/watertap/unit_models/translators/tests/test_translator_adm1_asm2d.py b/watertap/unit_models/translators/tests/test_translator_adm1_asm2d.py
index cf71b41b1e..5e674a2e6b 100644
--- a/watertap/unit_models/translators/tests/test_translator_adm1_asm2d.py
+++ b/watertap/unit_models/translators/tests/test_translator_adm1_asm2d.py
@@ -604,6 +604,4 @@ def test_example_case_scaler(self):
         # Check condition number to confirm scaling
         sm = TransformationFactory("core.scale_model").create_using(m, rename=False)
         jac, _ = get_jacobian(sm, scaled=False)
-        assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            7.180968e3, rel=1e-3
-        )
+        assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(712.33, rel=1e-3)
diff --git a/watertap/unit_models/translators/tests/test_translator_asm1_adm1.py b/watertap/unit_models/translators/tests/test_translator_asm1_adm1.py
index 085f46f01d..67aab650ef 100644
--- a/watertap/unit_models/translators/tests/test_translator_asm1_adm1.py
+++ b/watertap/unit_models/translators/tests/test_translator_asm1_adm1.py
@@ -700,5 +700,5 @@ def test_example_case_scaler(self):
         sm = TransformationFactory("core.scale_model").create_using(m, rename=False)
         jac, _ = get_jacobian(sm, scaled=False)
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            2.20621e5, rel=1e-3
+            6.11581e5, rel=1e-3
         )
diff --git a/watertap/unit_models/translators/tests/test_translator_asm2d_adm1.py b/watertap/unit_models/translators/tests/test_translator_asm2d_adm1.py
index aa69c95553..1e9b50c5b7 100644
--- a/watertap/unit_models/translators/tests/test_translator_asm2d_adm1.py
+++ b/watertap/unit_models/translators/tests/test_translator_asm2d_adm1.py
@@ -728,9 +728,7 @@ def test_constraint_scaling_routine(self, model):
 
         scaler.constraint_scaling_routine(model.fs.unit)
 
-        sfx_out = model.fs.unit.properties_out[0].scaling_factor
-        assert isinstance(sfx_out, Suffix)
-        assert len(sfx_out) == 0
+        assert not hasattr(model.fs.unit.properties_out[0], "scaling_factor")
 
     @pytest.mark.component
     def test_scale_model(self, model):
@@ -920,5 +918,5 @@ def test_example_case_scaler(self):
         sm = TransformationFactory("core.scale_model").create_using(m, rename=False)
         jac, _ = get_jacobian(sm, scaled=False)
         assert (jacobian_cond(jac=jac, scaled=False)) == pytest.approx(
-            5.640234319533e3, rel=1e-3
+            3969.50, rel=1e-3
         )