diff --git a/.github/workflows/checks.yml b/.github/workflows/checks.yml
index c3f5cecc24..dff736ef8c 100644
--- a/.github/workflows/checks.yml
+++ b/.github/workflows/checks.yml
@@ -258,6 +258,10 @@ jobs:
         echo '::group::Output of "idaes get-extensions" command'
         idaes get-extensions --verbose
         echo '::endgroup::'
+    - name: Exclude notebooks that cause errors on Windows
+      if: startswith(matrix.os, 'win')
+      run: |
+        rm tutorials/nawi_spring_meeting2023.ipynb
     - name: Run pytest with nbmake
       run:
         pytest --nbmake **/*.ipynb
diff --git a/setup.py b/setup.py
index 36c6bec173..26b63673e2 100644
--- a/setup.py
+++ b/setup.py
@@ -21,14 +21,14 @@
 long_description = (cwd / "README.md").read_text()
 
 SPECIAL_DEPENDENCIES_FOR_RELEASE = [
-    "idaes-pse==2.0.*",  # from PyPI
+    "idaes-pse==2.1.*",  # from PyPI
 ]
 
 SPECIAL_DEPENDENCIES_FOR_PRERELEASE = [
     # update with a tag from the nawi-hub/idaes-pse
     # when a version of IDAES newer than the latest stable release from PyPI
     # will become needed for the watertap development
-    "idaes-pse==2.0.*",
+    "idaes-pse==2.1.0",
 ]
 
 # Arguments marked as "Required" below must be included for upload to PyPI.
@@ -80,7 +80,7 @@
         # primary requirements for unit and property models
         # maintainers: switch to SPECIAL_DEPENDENCIES_FOR_RELEASE when cutting a release of watertap
         *SPECIAL_DEPENDENCIES_FOR_PRERELEASE,
-        "pyomo>=6.2,<6.6",  # (also needed for units in electrolyte database (edb))
+        "pyomo>=6.6.1",  # (also needed for units in electrolyte database (edb))
         # the following requirements are for the electrolyte database (edb)
         "pymongo>3",  # database interface
         "fastjsonschema",  # schema validation
@@ -103,6 +103,7 @@
             "mongomock",
             "pandas",
             "nbmake",
+            "nbconvert",
         ],
         "dev": [
             "nbsphinx",  # jupyter notebook support for sphinx
diff --git a/watertap/core/tests/test_zero_order_properties.py b/watertap/core/tests/test_zero_order_properties.py
index aac1f26acc..aa146b2bac 100644
--- a/watertap/core/tests/test_zero_order_properties.py
+++ b/watertap/core/tests/test_zero_order_properties.py
@@ -210,9 +210,9 @@ def test_CV_integration(model):
 
     model.fs.cv.add_geometry()
 
-    model.fs.cv.add_state_blocks(has_phase_equilibrium=True)
+    model.fs.cv.add_state_blocks(has_phase_equilibrium=False)
 
-    model.fs.cv.add_material_balances(has_phase_equilibrium=True)
+    model.fs.cv.add_material_balances(has_phase_equilibrium=False)
 
     # No energy or momentum balances, as these are not supported.
 
diff --git a/watertap/costing/watertap_costing_package.py b/watertap/costing/watertap_costing_package.py
index 47b6d0d12a..782e4e792e 100644
--- a/watertap/costing/watertap_costing_package.py
+++ b/watertap/costing/watertap_costing_package.py
@@ -274,13 +274,27 @@ def build_process_costs(self):
             == self.factor_maintenance_labor_chemical * self.total_capital_cost
         )
 
-        self.total_operating_cost_constraint = pyo.Constraint(
-            expr=self.total_operating_cost
-            == self.maintenance_labor_chemical_operating_cost
-            + self.aggregate_fixed_operating_cost
-            + self.aggregate_variable_operating_cost
-            + sum(self.aggregate_flow_costs.values()) * self.utilization_factor
-        )
+        if (
+            pyo.units.get_units(sum(self.aggregate_flow_costs.values()))
+        ) == pyo.units.dimensionless:
+            self.total_operating_cost_constraint = pyo.Constraint(
+                expr=self.total_operating_cost
+                == self.maintenance_labor_chemical_operating_cost
+                + self.aggregate_fixed_operating_cost
+                + self.aggregate_variable_operating_cost
+                + sum(self.aggregate_flow_costs.values())
+                * self.base_currency
+                / self.base_period
+                * self.utilization_factor
+            )
+        else:
+            self.total_operating_cost_constraint = pyo.Constraint(
+                expr=self.total_operating_cost
+                == self.maintenance_labor_chemical_operating_cost
+                + self.aggregate_fixed_operating_cost
+                + self.aggregate_variable_operating_cost
+                + sum(self.aggregate_flow_costs.values()) * self.utilization_factor
+            )
 
     def initialize_build(self):
         calculate_variable_from_constraint(
diff --git a/watertap/examples/chemistry/tests/test_pH_dependent_solubility.py b/watertap/examples/chemistry/tests/test_pH_dependent_solubility.py
index ccce31b334..2bf3b9f9d4 100644
--- a/watertap/examples/chemistry/tests/test_pH_dependent_solubility.py
+++ b/watertap/examples/chemistry/tests/test_pH_dependent_solubility.py
@@ -482,6 +482,7 @@ def test_case_1_no_dissolution():
     )
 
 
+@pytest.mark.requires_idaes_solver
 @pytest.mark.component
 def test_case_1_high_dissolution():
     model = run_case1(
@@ -494,6 +495,7 @@ def test_case_1_high_dissolution():
     )
 
 
+@pytest.mark.requires_idaes_solver
 @pytest.mark.component
 def test_case_1_mid_dissolution():
     model = run_case1(
@@ -506,6 +508,7 @@ def test_case_1_mid_dissolution():
     )
 
 
+@pytest.mark.requires_idaes_solver
 @pytest.mark.component
 def test_case_1_low_dissolution():
     model = run_case1(
diff --git a/watertap/examples/flowsheets/case_studies/activated_sludge/ASM2D_flowsheet.py b/watertap/examples/flowsheets/case_studies/activated_sludge/ASM2D_flowsheet.py
index f89cfd8621..21ba1242d6 100644
--- a/watertap/examples/flowsheets/case_studies/activated_sludge/ASM2D_flowsheet.py
+++ b/watertap/examples/flowsheets/case_studies/activated_sludge/ASM2D_flowsheet.py
@@ -282,13 +282,11 @@ def scale_variables(m):
             if "temperature" in var.name:
                 iscale.set_scaling_factor(var, 1e-1)
             if "pressure" in var.name:
-                iscale.set_scaling_factor(var, 1e-3)
+                iscale.set_scaling_factor(var, 1e-4)
             if "enth_mol" in var.name:
                 iscale.set_scaling_factor(var, 1e-3)
             if "alkalinity" in var.name:
                 iscale.set_scaling_factor(var, 1e3)
-            if "conc_mass_comp" in var.name:
-                iscale.set_scaling_factor(var, 1e2)
 
             if "conc_mass_comp[S_O2]" in var.name:
                 iscale.set_scaling_factor(var, 1e3)
@@ -398,7 +396,7 @@ def function(unit):
 
     seq.run(m, function)
 
-    solver = get_solver(options={"bound_push": 1e-8})
+    solver = get_solver()
     results = solver.solve(m, tee=False)
     check_solve(results, checkpoint="closing recycle", logger=_log, fail_flag=True)
 
diff --git a/watertap/examples/flowsheets/case_studies/activated_sludge/tests/test_asm2d_flowsheet.py b/watertap/examples/flowsheets/case_studies/activated_sludge/tests/test_asm2d_flowsheet.py
index 73a07f02be..6265e83b44 100644
--- a/watertap/examples/flowsheets/case_studies/activated_sludge/tests/test_asm2d_flowsheet.py
+++ b/watertap/examples/flowsheets/case_studies/activated_sludge/tests/test_asm2d_flowsheet.py
@@ -57,10 +57,10 @@ def test_results(self, model):
         assert value(model.fs.Treated.temperature[0]) == pytest.approx(298.15, rel=1e-4)
         assert value(model.fs.Treated.pressure[0]) == pytest.approx(101325, rel=1e-4)
         assert value(model.fs.Treated.conc_mass_comp[0, "S_A"]) == pytest.approx(
-            5.3211e-5, rel=1e-4
+            5.4560e-5, rel=1e-4
         )
         assert value(model.fs.Treated.conc_mass_comp[0, "S_F"]) == pytest.approx(
-            2.5618e-3, rel=1e-2
+            2.5736e-3, rel=1e-2
         )
         assert value(model.fs.Treated.conc_mass_comp[0, "S_I"]) == pytest.approx(
             30e-3, rel=1e-4
@@ -69,25 +69,25 @@ def test_results(self, model):
             13.837e-3, abs=1e-4
         )
         assert value(model.fs.Treated.conc_mass_comp[0, "S_NH4"]) == pytest.approx(
-            9.0286e-6, rel=1e-4
+            8.98349e-6, rel=1e-4
         )
         assert value(model.fs.Treated.conc_mass_comp[0, "S_NO3"]) == pytest.approx(
-            4.4784e-3, abs=1e-4
+            4.3924e-3, abs=1e-4
         )
         assert value(model.fs.Treated.conc_mass_comp[0, "S_O2"]) == pytest.approx(
-            7.8653e-3, abs=1e-4
+            7.8744e-3, abs=1e-4
         )
         assert value(model.fs.Treated.conc_mass_comp[0, "S_PO4"]) == pytest.approx(
-            2.5913e-3, rel=1e-4
+            2.7958e-3, rel=1e-4
         )
         assert value(model.fs.Treated.conc_mass_comp[0, "X_AUT"]) == pytest.approx(
-            1.3278e-3, abs=1e-4
+            1.3277e-3, abs=1e-4
         )
         assert value(model.fs.Treated.conc_mass_comp[0, "X_H"]) == pytest.approx(
-            36.313e-3, rel=1e-4
+            36.331e-3, rel=1e-4
         )
         assert value(model.fs.Treated.conc_mass_comp[0, "X_I"]) == pytest.approx(
-            21.178e-3, rel=1e-4
+            21.1808e-3, rel=1e-4
         )
         assert value(model.fs.Treated.conc_mass_comp[0, "X_MeOH"]) == pytest.approx(
             0, abs=1e-4
@@ -96,32 +96,32 @@ def test_results(self, model):
             0, abs=1e-4
         )
         assert value(model.fs.Treated.conc_mass_comp[0, "X_PAO"]) == pytest.approx(
-            2.21218e-3, abs=1e-4
+            2.111e-3, abs=1e-4
         )
         assert value(model.fs.Treated.conc_mass_comp[0, "X_PHA"]) == pytest.approx(
-            2.1218e-9, abs=1e-4
+            8.981e-5, abs=1e-4
         )
         assert value(model.fs.Treated.conc_mass_comp[0, "X_PP"]) == pytest.approx(
-            0.7773e-3, abs=1e-4
+            0.66913e-3, abs=1e-4
         )
         assert value(model.fs.Treated.conc_mass_comp[0, "X_S"]) == pytest.approx(
-            1.1586e-3, rel=1e-4
+            1.15897e-3, rel=1e-4
         )
         assert value(model.fs.Treated.conc_mass_comp[0, "X_TSS"]) == pytest.approx(
-            76.828e-3, rel=1e-4
+            76.489e-3, rel=1e-4
         )
         assert value(model.fs.Treated.alkalinity[0]) == pytest.approx(
-            7.2343e-3, rel=1e-4
+            7.243e-3, rel=1e-4
         )
         # Sludge stream
         assert value(model.fs.Sludge.flow_vol[0]) == pytest.approx(4.2808e-3, rel=1e-4)
         assert value(model.fs.Sludge.temperature[0]) == pytest.approx(298.15, rel=1e-4)
         assert value(model.fs.Sludge.pressure[0]) == pytest.approx(101325, rel=1e-4)
         assert value(model.fs.Sludge.conc_mass_comp[0, "S_A"]) == pytest.approx(
-            5.3211e-5, rel=1e-4
+            5.4560e-5, rel=1e-4
         )
         assert value(model.fs.Sludge.conc_mass_comp[0, "S_F"]) == pytest.approx(
-            2.5618e-3, rel=1e-2
+            2.5736e-3, rel=1e-2
         )
         assert value(model.fs.Sludge.conc_mass_comp[0, "S_I"]) == pytest.approx(
             30e-3, rel=1e-4
@@ -130,25 +130,25 @@ def test_results(self, model):
             13.837e-3, abs=1e-4
         )
         assert value(model.fs.Sludge.conc_mass_comp[0, "S_NH4"]) == pytest.approx(
-            9.0286e-6, rel=1e-4
+            8.9834e-6, rel=1e-4
         )
         assert value(model.fs.Sludge.conc_mass_comp[0, "S_NO3"]) == pytest.approx(
-            4.4784e-3, abs=1e-4
+            4.3924e-3, abs=1e-4
         )
         assert value(model.fs.Sludge.conc_mass_comp[0, "S_O2"]) == pytest.approx(
-            7.8653e-3, abs=1e-4
+            7.8744e-3, abs=1e-4
         )
         assert value(model.fs.Sludge.conc_mass_comp[0, "S_PO4"]) == pytest.approx(
-            2.5913e-3, rel=1e-4
+            2.7958e-3, rel=1e-4
         )
         assert value(model.fs.Sludge.conc_mass_comp[0, "X_AUT"]) == pytest.approx(
-            58.680e-3, abs=1e-4
+            58.6744e-3, abs=1e-4
         )
         assert value(model.fs.Sludge.conc_mass_comp[0, "X_H"]) == pytest.approx(
-            1.6193, rel=1e-4
+            1.6200, rel=1e-4
         )
         assert value(model.fs.Sludge.conc_mass_comp[0, "X_I"]) == pytest.approx(
-            953.56e-3, rel=1e-4
+            953.686e-3, rel=1e-4
         )
         assert value(model.fs.Sludge.conc_mass_comp[0, "X_MeOH"]) == pytest.approx(
             0, abs=1e-4
@@ -157,20 +157,20 @@ def test_results(self, model):
             0, abs=1e-4
         )
         assert value(model.fs.Sludge.conc_mass_comp[0, "X_PAO"]) == pytest.approx(
-            94.272e-3, abs=1e-4
+            93.795e-3, abs=1e-4
         )
         assert value(model.fs.Sludge.conc_mass_comp[0, "X_PHA"]) == pytest.approx(
-            3.6979e-3, abs=1e-4
+            4.0913e-3, abs=1e-4
         )
         assert value(model.fs.Sludge.conc_mass_comp[0, "X_PP"]) == pytest.approx(
-            35.460e-3, abs=1e-4
+            30.522e-3, abs=1e-4
         )
         assert value(model.fs.Sludge.conc_mass_comp[0, "X_S"]) == pytest.approx(
-            50.980e-3, rel=1e-4
+            50.995e-3, rel=1e-4
         )
         assert value(model.fs.Sludge.conc_mass_comp[0, "X_TSS"]) == pytest.approx(
-            3.423, rel=1e-4
+            3.4078, rel=1e-4
         )
         assert value(model.fs.Sludge.alkalinity[0]) == pytest.approx(
-            7.2343e-3, rel=1e-4
+            7.2430e-3, rel=1e-4
         )
diff --git a/watertap/examples/flowsheets/case_studies/electroNP/electroNP_flowsheet.py b/watertap/examples/flowsheets/case_studies/electroNP/electroNP_flowsheet.py
index ed2a4f74f8..fc4d0771c5 100644
--- a/watertap/examples/flowsheets/case_studies/electroNP/electroNP_flowsheet.py
+++ b/watertap/examples/flowsheets/case_studies/electroNP/electroNP_flowsheet.py
@@ -213,6 +213,32 @@ def build_flowsheet():
     iscale.calculate_scaling_factors(m)
 
     iscale.set_scaling_factor(m.fs.electroNP.properties_byproduct[0.0].flow_vol, 1e7)
+    iscale.set_scaling_factor(m.fs.AD.vapor_phase[0].pressure_sat, 1e-3)
+
+    skip_constraint = [
+        "S_O2",
+        "S_N2",
+        "X_AUT",
+        "X_PHA",
+        "S_I",
+        "S_NO3",
+        "S_Mg",
+        "X_MeP",
+        "X_MeOH",
+        "X_PAO",
+        "X_TSS",
+        "S_F",
+        "S_K",
+        "S_A",
+        "X_I",
+        "X_H",
+        "X_PP",
+        "X_S",
+    ]
+
+    for i in skip_constraint:
+        m.fs.electroNP.solute_removal_equation[0.0, "Liq", i].deactivate()
+        m.fs.electroNP.solute_treated_equation[0.0, "Liq", i].deactivate()
 
     m.fs.AD.initialize(outlvl=idaeslog.INFO_HIGH)
     propagate_state(m.fs.stream_adm1_translator)
diff --git a/watertap/examples/flowsheets/case_studies/wastewater_resource_recovery/amo_1690/amo_1690.py b/watertap/examples/flowsheets/case_studies/wastewater_resource_recovery/amo_1690/amo_1690.py
index d90e29e413..ee42e88653 100644
--- a/watertap/examples/flowsheets/case_studies/wastewater_resource_recovery/amo_1690/amo_1690.py
+++ b/watertap/examples/flowsheets/case_studies/wastewater_resource_recovery/amo_1690/amo_1690.py
@@ -142,6 +142,8 @@ def set_operating_conditions(m):
     m.fs.feed.conc_mass_comp[0, "acetic_acid"].fix(conc_mass_aa)
     m.fs.feed.conc_mass_comp[0, "ammonium_as_nitrogen"].fix(conc_mass_nh4)
 
+    m.fs.feed.initialize()
+
     solve(m.fs.feed, checkpoint="solve feed block")
 
     m.fs.cmf.load_parameters_from_database(use_default_removal=True)
diff --git a/watertap/examples/flowsheets/lsrro/lsrro.py b/watertap/examples/flowsheets/lsrro/lsrro.py
index eca63565e8..79865f2660 100644
--- a/watertap/examples/flowsheets/lsrro/lsrro.py
+++ b/watertap/examples/flowsheets/lsrro/lsrro.py
@@ -406,7 +406,7 @@ def stage_recovery_mass_H2O(fs, stage):
             + (
                 m.fs.costing.booster_pump_capex_lcow
                 if number_of_stages > 1
-                else 0 * m.fs.costing.base_currency
+                else 0 * m.fs.costing.base_currency / pyunits.m**3
             )
             + m.fs.costing.erd_capex_lcow
         )
diff --git a/watertap/examples/flowsheets/nf_dspmde/tests/test_nf_with_bypass.py b/watertap/examples/flowsheets/nf_dspmde/tests/test_nf_with_bypass.py
index 979f556c0a..b731082998 100644
--- a/watertap/examples/flowsheets/nf_dspmde/tests/test_nf_with_bypass.py
+++ b/watertap/examples/flowsheets/nf_dspmde/tests/test_nf_with_bypass.py
@@ -9,12 +9,17 @@
 # information, respectively. These files are also available online at the URL
 # "https://github.com/watertap-org/watertap/"
 #################################################################################
+import sys
 
 import pytest
 from pyomo.environ import value
 from watertap.examples.flowsheets.nf_dspmde.nf_with_bypass import main
 
 
+@pytest.mark.skipif(
+    sys.platform.startswith("win"),
+    reason="Known issue on Windows, see watertap-org/watertap#1072",
+)
 @pytest.mark.requires_idaes_solver
 @pytest.mark.component
 def test_main():
diff --git a/watertap/property_models/NDMA_prop_pack.py b/watertap/property_models/NDMA_prop_pack.py
index 509aa922ba..0cf5950398 100644
--- a/watertap/property_models/NDMA_prop_pack.py
+++ b/watertap/property_models/NDMA_prop_pack.py
@@ -159,6 +159,21 @@ class _NDMAStateBlock(StateBlock):
     whole, rather than individual elements of indexed Property Blocks.
     """
 
+    def fix_initialization_states(self):
+        """
+        Fixes state variables for state blocks.
+
+        Returns:
+            None
+        """
+        # Fix state variables
+        fix_state_vars(self)
+
+        # Constraint on water concentration at outlet - unfix in these cases
+        for b in self.values():
+            if b.config.defined_state is False:
+                b.conc_mol_comp["H2O"].unfix()
+
     def initialize(
         self,
         state_args=None,
diff --git a/watertap/property_models/NaCl_T_dep_prop_pack.py b/watertap/property_models/NaCl_T_dep_prop_pack.py
index ed90915aaa..e9688dd024 100644
--- a/watertap/property_models/NaCl_T_dep_prop_pack.py
+++ b/watertap/property_models/NaCl_T_dep_prop_pack.py
@@ -726,6 +726,21 @@ class _NaClStateBlock(StateBlock):
     whole, rather than individual elements of indexed Property Blocks.
     """
 
+    def fix_initialization_states(self):
+        """
+        Fixes state variables for state blocks.
+
+        Returns:
+            None
+        """
+        # Fix state variables
+        fix_state_vars(self)
+
+        # Constraint on water concentration at outlet - unfix in these cases
+        for b in self.values():
+            if b.config.defined_state is False:
+                b.conc_mol_comp["H2O"].unfix()
+
     def initialize(
         self,
         state_args=None,
diff --git a/watertap/property_models/NaCl_prop_pack.py b/watertap/property_models/NaCl_prop_pack.py
index d75972fad7..3fc14c24d6 100644
--- a/watertap/property_models/NaCl_prop_pack.py
+++ b/watertap/property_models/NaCl_prop_pack.py
@@ -214,6 +214,21 @@ class _NaClStateBlock(StateBlock):
     whole, rather than individual elements of indexed Property Blocks.
     """
 
+    def fix_initialization_states(self):
+        """
+        Fixes state variables for state blocks.
+
+        Returns:
+            None
+        """
+        # Fix state variables
+        fix_state_vars(self)
+
+        # Constraint on water concentration at outlet - unfix in these cases
+        for b in self.values():
+            if b.config.defined_state is False:
+                b.conc_mol_comp["H2O"].unfix()
+
     def initialize(
         self,
         state_args=None,
diff --git a/watertap/property_models/activated_sludge/tests/test_asm1_integration.py b/watertap/property_models/activated_sludge/tests/test_asm1_integration.py
index c881613f32..9448a77a8d 100644
--- a/watertap/property_models/activated_sludge/tests/test_asm1_integration.py
+++ b/watertap/property_models/activated_sludge/tests/test_asm1_integration.py
@@ -12,7 +12,7 @@
 """
 Tests for ASM1 property and reaction packages.
 
-Verified againts results from:
+Verified against results from:
 
 [1] J. Alex, L. Benedetti, J. Copp, K.V. Gernaey, U. Jeppsson, I. Nopens, M.N. Pons,
 J.P. Steyer and P. Vanrolleghem, "Benchmark Simulation Model no. 1 (BSM1)", 2018
diff --git a/watertap/property_models/activated_sludge/tests/test_asm1_reaction.py b/watertap/property_models/activated_sludge/tests/test_asm1_reaction.py
index 4b3be99cd1..6b476100b0 100644
--- a/watertap/property_models/activated_sludge/tests/test_asm1_reaction.py
+++ b/watertap/property_models/activated_sludge/tests/test_asm1_reaction.py
@@ -11,6 +11,7 @@
 #################################################################################
 """
 Tests for ASM1 reaction package.
+
 Authors: Andrew Lee
 """
 import pytest
diff --git a/watertap/property_models/cryst_prop_pack.py b/watertap/property_models/cryst_prop_pack.py
index f66ed5459c..a8a22744a5 100644
--- a/watertap/property_models/cryst_prop_pack.py
+++ b/watertap/property_models/cryst_prop_pack.py
@@ -874,6 +874,21 @@ class _NaClStateBlock(StateBlock):
     whole, rather than individual elements of indexed Property Blocks.
     """
 
+    def fix_initialization_states(self):
+        """
+        Fixes state variables for state blocks.
+
+        Returns:
+            None
+        """
+        # Fix state variables
+        fix_state_vars(self)
+
+        # Constraint on water concentration at outlet - unfix in these cases
+        for b in self.values():
+            if b.config.defined_state is False:
+                b.conc_mol_comp["H2O"].unfix()
+
     def initialize(
         self,
         state_args=None,
diff --git a/watertap/property_models/seawater_prop_pack.py b/watertap/property_models/seawater_prop_pack.py
index 02f8a11512..493d99215e 100644
--- a/watertap/property_models/seawater_prop_pack.py
+++ b/watertap/property_models/seawater_prop_pack.py
@@ -784,6 +784,21 @@ class _SeawaterStateBlock(StateBlock):
     whole, rather than individual elements of indexed Property Blocks.
     """
 
+    def fix_initialization_states(self):
+        """
+        Fixes state variables for state blocks.
+
+        Returns:
+            None
+        """
+        # Fix state variables
+        fix_state_vars(self)
+
+        # Constraint on water concentration at outlet - unfix in these cases
+        for b in self.values():
+            if b.config.defined_state is False:
+                b.conc_mol_comp["H2O"].unfix()
+
     def initialize(
         self,
         state_args=None,
diff --git a/watertap/property_models/water_prop_pack.py b/watertap/property_models/water_prop_pack.py
index cd2f9bba56..ad622ad25f 100644
--- a/watertap/property_models/water_prop_pack.py
+++ b/watertap/property_models/water_prop_pack.py
@@ -398,6 +398,21 @@ class _WaterStateBlock(StateBlock):
     whole, rather than individual elements of indexed Property Blocks.
     """
 
+    def fix_initialization_states(self):
+        """
+        Fixes state variables for state blocks.
+
+        Returns:
+            None
+        """
+        # Fix state variables
+        fix_state_vars(self)
+
+        # Constraint on water concentration at outlet - unfix in these cases
+        for b in self.values():
+            if b.config.defined_state is False:
+                b.conc_mol_comp["H2O"].unfix()
+
     def initialize(
         self,
         state_args=None,
diff --git a/watertap/unit_models/tests/test_crystallizer.py b/watertap/unit_models/tests/test_crystallizer.py
index 005446c896..c304dd45cd 100644
--- a/watertap/unit_models/tests/test_crystallizer.py
+++ b/watertap/unit_models/tests/test_crystallizer.py
@@ -322,6 +322,7 @@ def test_calculate_scaling(self, Crystallizer_frame):
         for _ in badly_scaled_var_generator(m):
             assert False
 
+    @pytest.mark.requires_idaes_solver
     @pytest.mark.component
     def test_initialize(self, Crystallizer_frame):
         # Add costing function, then initialize
@@ -342,6 +343,7 @@ def test_initialize(self, Crystallizer_frame):
     #     badly_scaled_var_lst = list(badly_scaled_var_generator(m))
     #     assert badly_scaled_var_lst == []
 
+    @pytest.mark.requires_idaes_solver
     @pytest.mark.component
     def test_solve(self, Crystallizer_frame):
         m = Crystallizer_frame
@@ -351,6 +353,7 @@ def test_solve(self, Crystallizer_frame):
         assert results.solver.termination_condition == TerminationCondition.optimal
         assert results.solver.status == SolverStatus.ok
 
+    @pytest.mark.requires_idaes_solver
     @pytest.mark.component
     def test_conservation(self, Crystallizer_frame):
         m = Crystallizer_frame
@@ -411,6 +414,7 @@ def test_conservation(self, Crystallizer_frame):
             <= 1e-2
         )
 
+    @pytest.mark.requires_idaes_solver
     @pytest.mark.component
     def test_solution(self, Crystallizer_frame):
         m = Crystallizer_frame
@@ -454,6 +458,7 @@ def test_solution(self, Crystallizer_frame):
             m.fs.costing.aggregate_capital_cost
         )
 
+    @pytest.mark.requires_idaes_solver
     @pytest.mark.component
     def test_solution2_capcosting_by_mass(self, Crystallizer_frame):
         m = Crystallizer_frame