diff --git a/watertap/costing/units/crystallizer.py b/watertap/costing/units/crystallizer.py
index 42971184d7..f8f66013a6 100644
--- a/watertap/costing/units/crystallizer.py
+++ b/watertap/costing/units/crystallizer.py
@@ -86,8 +86,15 @@ def build_crystallizer_cost_param_block(blk):
         doc="Steam cost, Panagopoulos (2019)",
     )
 
+    blk.NaCl_recovery_value = pyo.Var(
+        initialize=0,
+        units=pyo.units.USD_2018 / pyo.units.kg,
+        doc="Unit recovery value of NaCl",
+    )
+
     costing = blk.parent_block()
     costing.add_defined_flow("steam", blk.steam_cost)
+    costing.add_defined_flow("NaCl", blk.NaCl_recovery_value)
 
 
 def cost_crystallizer(blk, cost_type=CrystallizerCostType.default):
@@ -135,6 +142,11 @@ def _cost_crystallizer_flows(blk):
         "steam",
     )
 
+    blk.costing_package.cost_flow(
+        blk.unit_model.solids.flow_mass_phase_comp[0, "Sol", "NaCl"],
+        "NaCl",
+    )
+
 
 @register_costing_parameter_block(
     build_rule=build_crystallizer_cost_param_block,
diff --git a/watertap/costing/watertap_costing_package.py b/watertap/costing/watertap_costing_package.py
index 6ed3bd742c..7d66b3629e 100644
--- a/watertap/costing/watertap_costing_package.py
+++ b/watertap/costing/watertap_costing_package.py
@@ -252,7 +252,7 @@ def build_process_costs(self):
         )
         self.total_operating_cost = pyo.Var(
             initialize=1e3,
-            domain=pyo.NonNegativeReals,
+            domain=pyo.Reals,
             doc="Total operating cost",
             units=self.base_currency / self.base_period,
         )
diff --git a/watertap/examples/flowsheets/full_treatment_train/analysis/tests/test_flowsheets.py b/watertap/examples/flowsheets/full_treatment_train/analysis/tests/test_flowsheets.py
index 25bab6ab4b..20282416df 100644
--- a/watertap/examples/flowsheets/full_treatment_train/analysis/tests/test_flowsheets.py
+++ b/watertap/examples/flowsheets/full_treatment_train/analysis/tests/test_flowsheets.py
@@ -30,7 +30,7 @@ def test_flowsheet_NF():
     ) == pytest.approx(0.7667, rel=1e-3)
     assert value(
         m.fs.tb_pretrt_to_desal.properties_in[0].flow_mass_phase_comp["Liq", "Ca"]
-    ) == pytest.approx(1.15808e-4, rel=1e-3)
+    ) == pytest.approx(1.1636e-4, rel=1e-3)
 
 
 @pytest.mark.component
diff --git a/watertap/examples/flowsheets/full_treatment_train/flowsheet_components/chemistry/tests/test_posttreatment.py b/watertap/examples/flowsheets/full_treatment_train/flowsheet_components/chemistry/tests/test_posttreatment.py
index 1d85d48a68..954eaa3426 100644
--- a/watertap/examples/flowsheets/full_treatment_train/flowsheet_components/chemistry/tests/test_posttreatment.py
+++ b/watertap/examples/flowsheets/full_treatment_train/flowsheet_components/chemistry/tests/test_posttreatment.py
@@ -72,24 +72,24 @@ def test_ideal_naocl_chlorination():
 def test_ideal_naocl_chlorination_full_block():
     model = run_chlorination_block_example(fix_free_chlorine=True)
     assert model.fs.ideal_naocl_mixer_unit.dosing_rate.value == pytest.approx(
-        1.7296113311683092e-09, rel=1e-3
+        9.505698559578499e-07, rel=1e-3
     )
     assert model.fs.ideal_naocl_mixer_unit.outlet.flow_mol[0].value == pytest.approx(
         25.000025535888078, rel=1e-3
     )
     assert model.fs.ideal_naocl_mixer_unit.outlet.mole_frac_comp[
         0, "OCl_-"
-    ].value == pytest.approx(1.6123004572288052e-07, rel=1e-3)
+    ].value == pytest.approx(5.107133802822922e-07, rel=1e-3)
 
     assert model.fs.ideal_naocl_chlorination_unit.free_chlorine.value == pytest.approx(
         2, rel=1e-3
     )
     assert model.fs.ideal_naocl_chlorination_unit.outlet.mole_frac_comp[
         0, "OCl_-"
-    ].value == pytest.approx(5.0027242332010015e-08, rel=1e-3)
+    ].value == pytest.approx(4.5088044031726496e-07, rel=1e-3)
     assert model.fs.ideal_naocl_chlorination_unit.outlet.mole_frac_comp[
         0, "H_+"
-    ].value == pytest.approx(1.49011416785194e-10, rel=1e-3)
+    ].value == pytest.approx(5.4260427865375997e-11, rel=1e-3)
 
 
 @pytest.mark.component
diff --git a/watertap/examples/flowsheets/nf_dspmde/tests/test_nf.py b/watertap/examples/flowsheets/nf_dspmde/tests/test_nf.py
index 18fb6787a6..7cf5bb32ea 100644
--- a/watertap/examples/flowsheets/nf_dspmde/tests/test_nf.py
+++ b/watertap/examples/flowsheets/nf_dspmde/tests/test_nf.py
@@ -20,12 +20,12 @@
 def test_main():
     m = main()
     test_dict = {
-        "lcow": [m.fs.costing.LCOW, 0.16811587158493219],
-        "pressure": [m.fs.NF.pump.outlet.pressure[0] / 1e5, 6.56],
-        "area": [m.fs.NF.nfUnit.area, 285.6900547389303],
+        "lcow": [m.fs.costing.LCOW, 0.15058960529129017],
+        "pressure": [m.fs.NF.pump.outlet.pressure[0] / 1e5, 8.13],
+        "area": [m.fs.NF.nfUnit.area, 423.8956418211484],
         "recovery": [
             m.fs.NF.nfUnit.recovery_vol_phase[0.0, "Liq"] * 100,
-            73.47934090302432,
+            94.99999441324391,
         ],
     }
     for (model_result, testval) in test_dict.values():
diff --git a/watertap/unit_models/tests/test_crystallizer.py b/watertap/unit_models/tests/test_crystallizer.py
index c304dd45cd..fa4c72ef83 100644
--- a/watertap/unit_models/tests/test_crystallizer.py
+++ b/watertap/unit_models/tests/test_crystallizer.py
@@ -567,6 +567,9 @@ def test_solution2_operatingcost(self, Crystallizer_frame_2):
         assert pytest.approx(30666.67, rel=1e-3) == value(
             m.fs.costing.aggregate_flow_costs["steam"]
         )
+        assert pytest.approx(0, rel=1e-3) == value(
+            m.fs.costing.aggregate_flow_costs["NaCl"]
+        )
 
     @pytest.mark.component
     def test_solution2_operatingcost_steampressure(self, Crystallizer_frame_2):
@@ -589,3 +592,29 @@ def test_solution2_operatingcost_steampressure(self, Crystallizer_frame_2):
         assert pytest.approx(21451.91, rel=1e-3) == value(
             m.fs.costing.aggregate_flow_costs["steam"]
         )
+        assert pytest.approx(0, abs=1e-6) == value(
+            m.fs.costing.aggregate_flow_costs["NaCl"]
+        )
+
+    @pytest.mark.component
+    def test_solution2_operatingcost_NaCl_revenue(self, Crystallizer_frame_2):
+        m = Crystallizer_frame_2
+        m.fs.costing.crystallizer.steam_pressure.fix(3)
+        m.fs.costing.crystallizer.NaCl_recovery_value.fix(-0.07)
+
+        results = solver.solve(m)
+
+        # Check for optimal solution
+        assert results.solver.termination_condition == TerminationCondition.optimal
+        assert results.solver.status == SolverStatus.ok
+
+        # Operating cost validation
+        assert pytest.approx(835.41, rel=1e-3) == value(
+            m.fs.costing.aggregate_flow_costs["electricity"]
+        )
+        assert pytest.approx(30666.67, rel=1e-3) == value(
+            m.fs.costing.aggregate_flow_costs["steam"]
+        )
+        assert pytest.approx(-187858.2, rel=1e-3) == value(
+            m.fs.costing.aggregate_flow_costs["NaCl"]
+        )