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 51ea6ec635..3923203a16 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
@@ -29,7 +29,6 @@
     UnitModelBlockData,
 )
 from idaes.models.unit_models import (
-    CSTR,
     Feed,
     Separator,
     Product,
@@ -68,6 +67,7 @@
 from idaes.models.unit_models.mixer import MomentumMixingType
 from watertap.unit_models.translators.translator_asm2d_adm1 import Translator_ASM2d_ADM1
 from watertap.unit_models.anaerobic_digester import AD
+from watertap.unit_models.cstr import CSTR
 from watertap.unit_models.dewatering import (
     DewateringUnit,
     ActivatedSludgeModelType as dewater_type,
@@ -95,7 +95,7 @@
 
 def main(bio_P=False):
     m = build(bio_P=bio_P)
-    set_operating_conditions(m)
+    set_operating_conditions(m, bio_P=bio_P)
 
     for mx in m.fs.mixers:
         mx.pressure_equality_constraints[0.0, 2].deactivate()
@@ -411,7 +411,7 @@ def mass_transfer_R7(self, t):
     return m
 
 
-def set_operating_conditions(m):
+def set_operating_conditions(m, bio_P=False):
     # Feed Water Conditions
     print(f"DOF before feed: {degrees_of_freedom(m)}")
     m.fs.FeedWater.flow_vol.fix(20935.15 * pyo.units.m**3 / pyo.units.day)
@@ -537,6 +537,10 @@ def scale_variables(m):
             if "conc_mass_comp" in var.name:
                 iscale.set_scaling_factor(var, 1e1)
 
+    for unit in ("R1", "R2", "R3", "R4"):
+        block = getattr(m.fs, unit)
+        iscale.set_scaling_factor(block.hydraulic_retention_time, 1e-3)
+
     for unit in ("R1", "R2", "R3", "R4", "R5", "R6", "R7"):
         block = getattr(m.fs, unit)
         iscale.set_scaling_factor(
@@ -550,7 +554,18 @@ def scale_variables(m):
 
     iscale.set_scaling_factor(m.fs.AD.KH_co2, 1e1)
     iscale.set_scaling_factor(m.fs.AD.KH_ch4, 1e1)
-    iscale.set_scaling_factor(m.fs.AD.KH_h2, 1e1)
+    iscale.set_scaling_factor(m.fs.AD.KH_h2, 1e2)
+
+    if bio_P:
+        iscale.set_scaling_factor(m.fs.AD.liquid_phase.heat, 1e3)
+        iscale.constraint_scaling_transform(
+            m.fs.AD.liquid_phase.enthalpy_balances[0], 1e-6
+        )
+    else:
+        iscale.set_scaling_factor(m.fs.AD.liquid_phase.heat, 1e2)
+        iscale.constraint_scaling_transform(
+            m.fs.AD.liquid_phase.enthalpy_balances[0], 1e-3
+        )
 
     # Apply scaling
     scale_variables(m)
@@ -738,6 +753,11 @@ def add_costing(m):
         if isinstance(block, UnitModelBlockData) and hasattr(block, "costing"):
             iscale.set_scaling_factor(block.costing.capital_cost, 1e-5)
 
+    iscale.constraint_scaling_transform(m.fs.AD.costing.capital_cost_constraint, 1e-6)
+    iscale.constraint_scaling_transform(
+        m.fs.dewater.costing.capital_cost_constraint, 1e-6
+    )
+
 
 def display_costing(m):
     print("Levelized cost of water: %.2f $/m3" % pyo.value(m.fs.costing.LCOW))
@@ -952,7 +972,7 @@ def display_performance_metrics(m):
 
 
 if __name__ == "__main__":
-    m, results = main(bio_P=False)
+    m, results = main(bio_P=True)
 
     stream_table = create_stream_table_dataframe(
         {
diff --git a/watertap/flowsheets/full_water_resource_recovery_facility/BSM2_P_extension_ui.py b/watertap/flowsheets/full_water_resource_recovery_facility/BSM2_P_extension_ui.py
index 5f5b1e4064..527e064877 100644
--- a/watertap/flowsheets/full_water_resource_recovery_facility/BSM2_P_extension_ui.py
+++ b/watertap/flowsheets/full_water_resource_recovery_facility/BSM2_P_extension_ui.py
@@ -3873,7 +3873,7 @@ def build_flowsheet(build_options=None, **kwargs):
 
         m = build(bio_P=bioP)
 
-        set_operating_conditions(m)
+        set_operating_conditions(m, bio_P=bioP)
 
         for mx in m.fs.mixers:
             mx.pressure_equality_constraints[0.0, 2].deactivate()
@@ -3912,7 +3912,7 @@ def build_flowsheet(build_options=None, **kwargs):
     else:
         m = build(bio_P=False)
 
-        set_operating_conditions(m)
+        set_operating_conditions(m, bio_P=False)
 
         for mx in m.fs.mixers:
             mx.pressure_equality_constraints[0.0, 2].deactivate()
diff --git a/watertap/unit_models/anaerobic_digester.py b/watertap/unit_models/anaerobic_digester.py
index ec4c2f09eb..a297b6dcee 100644
--- a/watertap/unit_models/anaerobic_digester.py
+++ b/watertap/unit_models/anaerobic_digester.py
@@ -838,7 +838,6 @@ def calculate_scaling_factors(self):
 
         # TODO: improve this later; for now, this resolved some scaling issues for modified adm1 test file
         if "S_IP" in self.config.liquid_property_package.component_list:
-            iscale.set_scaling_factor(self.liquid_phase.heat, 1e-6)
             sf = iscale.get_scaling_factor(
                 self.liquid_phase.properties_out[0].conc_mass_comp["S_IP"],
                 default=1e-5,